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ABSTRACT 

The physical processes that heat the solar corona and accelerate the solar wind remain unknown after many 
years of study. Some have suggested that the wind is driven by waves and turbulence in open magnetic flux 
tubes, and others have suggested that plasma is injected into the open tubes by magnetic reconnection with 
closed loops. In order to test the latter idea, we developed Monte Carlo simulations of the photospheric "mag- 
netic carpet" and extrapolated the time-varying coronal field. These models were constructed for a range of 
different magnetic flux imbalance ratios. Completely balanced models represent quiet regions on the Sun and 
source regions of slow solar wind streams. Highly imbalanced models represent coronal holes and source 
regions of fast wind streams. The models agree with observed emergence rates, surface flux densities, and 
number distributions of magnetic elements. Despite having no imposed supergranular motions in the models, 
a realistic network of magnetic "funnels" appeared spontaneously. We computed the rate at which closed field 
lines open up (i.e., recycling times for open flux), and we estimated the energy flux released in reconnection 
events involving the opening up of closed flux tubes. For quiet regions and mixed-polarity coronal holes, these 
energy fluxes were found to be much lower than required to accelerate the solar wind. For the most imbalanced 
coronal holes, the energy fluxes may be large enough to power the solar wind, but the recycling times are far 
longer than the time it takes the solar wind to accelerate into the low corona. Thus, it is unlikely that either the 
slow or fast solar wind is driven by reconnection and loop-opening processes in the magnetic carpet. 
Subject headings: magnetic fields — magnetohydrodynamics (MHD) — plasmas — solar wind — Sun: corona 
— Sun: photosphere 



1. INTRODUCTION 

The magnetic field in the solar photosphere exists in a com- 
plex and continually evolving state that is driven by convec- 
tive motions under the surface. The dynamic interplay be- 
tween the magnetic field and the plasma has been called the 
Sun's "magnetic carpet" (Title and Schriiver 1998). There is 
a clear correlation between the topology and strength of the 
magnetic field and the energy deposition that is responsible 
for the hot (T > 10^ K) solar corona. We also know that the 
gas pressure associated with coronal heating is an important 
contributor to accelerating the supersonic solar wind (Parker 
[1958). Thus, it is natural to wonder to what extent the mag- 
netohydrodynamic (MHD) motions in the magnetic carpet are 
ultimately responsible for producing at least some of the solar 
wind's mass loss. 

Recently, two distinct classes of theoretical explanation 
have been proposed for the combined problem of coronal 
heating and solar wind acceleration. In the wave/turbulence- 
driven (WTD) models, convection jostles the open magnetic 
flux tubes that are rooted in the photosphere and produces 
waves that propagate into the corona. These waves (usually 
assumed to be Alfven waves) are proposed to partially 
reflect back down toward the Sun, develop into MHD 
turbule nce, and heat the plasma by th e ir gradual dissi- 
pation (iHollwed [T986t IVelli et al.1 [19911 IWanp & Sheeleyl 
199 It iMatthaeus et all Il999t ISuzuki & Inutsukal 12006 



occurs on small, supergranular scales (lAxfor d & McKenzie 
1991 iFisk etin [T99lTFisia I200I ISchwadron & McComas 
2003h . However, other models have been proposed in 



Cranmer et al.' '2007: 'Wang et al.r i2009l: IVerdini et al] I20T0 : 
Matsumoto & Shibata 2010). In the reconnection/loop- 
opening (RLO) class of models, it is assumed that closed, 
loop-like magnetic flux systems are the dominant source 
of mass and energy into the open-field regions. Some 
have suggested that RLO-type energy exchange primarily 
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which the "interchange reconnection" occurs in and between 
large- scale coronal s treamers further from t he solar surface 
(Eina udietal.lll999t ISuess &Nerneyll2004t lAntiochos et"aD 
i201Q) . 

The WTD idea of a flux tube that is open — and which stays 
open as the wind accelerates — is conceptually simpler than 
the idea of frequent changes in the flux tube topology. Be- 
cause of this simplicity, the WTD models have been subject 
to a greater degree of development and testing than the RLO 
models. In addition, we have a great deal of observational 
evidence that waves and turbulent motions are present ev- 
erywhere from the p hotosphere to the helio s phere (see, e.g., 
Tu & Marsch 1995; Bruno & Carbonell2005t Hans teen 2007] 
Aschwandenii2008i) . Thus, it is of interest to pursue the WTD 
idea to see how these waves affect the mean state of the 
plasma i n the absence o f any other sources of energy. For 
example. ICranmer et al.l (12007) and Cranmer (2009) showed 
that a set of WTD models that varied only the magnetic flux- 
tube expansion rate (and kept all other parameters fixed, in- 
cluding the wave fluxes at the lower boundary) can success- 
fully predict a wide range of measured properties of both fast 
and slow solar wind streams. 

RLO models need to be subjected to the same degree of 
development, testing, and refinement as the WTD models. 
This idea has a natural appeal since the op en flux tubes 
must be rooted in the vicinity of closed loops dDowdv et al.l 
Il986l) . In fact, multiple RLO-like reconnection events have 
been observed in coronal holes as "polar je ts" by instruments 
aboard SOHO Hinode, and STEREO (e .g.,|Wang.eLaD[1998t 
IShimojo et al.ll2007HNistic6 et al.ll2009l) . Reconnection at flie 
edges of coronal holes may be necessary to produce their ob- 
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served rigid rotation dLionello et al.ll2006l) . There are also ob- 
served correlations between the lengths of coronal loops, the 
electron temperature in the low cor ona, and the wind speed 
in interplanetary space (iGloeckler e t al. 2003) that are highly 
suggestive of a net transfer of magnetic energy from the loops 
to the open-field regions (see also Fi sk et al. 1999: Fisk 2003) . 

Testing the RLO idea using theoretical models is more diffi- 
cult than testing the WTD idea because of the complex multi- 
scale nature of the relevant magnetic fields. Many aspects 
of RLO-type processes cannot be simulated without resorting 
to fully three-dimensional and time-dependent models of the 
connection between the magnetic carpet and the solar wind. 
The goal of this paper is to begin constructing such mod- 
els in order to address several of the following unanswered 
questions about the RLO model. For example, how much of 
the magnetic energy that is liberated by reconnection goes 
into simply reconfiguring the closed fields, and how much 
goes into changing closed fields into open fields? Specifi- 
cally, what is the actual rate at which magnetic flux opens up 
from the magnetic carpet? Can the observed polar jets provide 
enough energy to drive a significant fraction of the solar wind? 
Lastly, how is the reconnection energy distributed into various 
forms (e.g., bulk kinetic energy, thermal energy, waves, or en- 
ergetic particles) that can each affect the accelerating wind in 
different ways? 

In this paper we present Monte Carlo models of the solar 
magnetic carpet that are used to determine the topology, tem- 
poral variability, and energy flux along field lines connected 
with the accelerating solar wind. Section|2]gives an overview 
of the motivations behind our choices of modeling technique. 
In Section [3] we describe the physical ingredients that went 
into the Monte Carlo models of the photospheric magnetic 
field. Section|4]then presents the results of these models and 
compares them with a range of observational diagnostics. In 
Section |5] we then describe how field lines were extrapolated 
from the photospheric lower boundary up into the corona, and 
we discuss the resulting time scales and energy fluxes that 
were derived for flux tubes relevant to RLO wind acceleration 
models. Finally, Section |6] concludes this paper with a brief 
summary of the major results, a discussion of some of the 
wider implications of this work, and suggestions for future 
improvements. 

2. MOTIVATIONS AND METHODS 

In this section we summarize the techniques that we chose 
to simulate the connections between the photospheric mag- 
netic field and the open flux tubes feeding the solar wind. 
It is also important to clarify how and why our assump- 
tions are consistent with the goal to quantify the impact of 
RLO physical processes. Our modeling was done in two 
steps. First, we simulated the photospheric magnetic carpet 
by means of a Monte Carlo ensemble of positive and negative 
monopole sources of magnetic flux. These sources are as- 
sumed to emerge from below (as bipolar ephemeral regions), 
move around on the surface, merge or cancel with their neigh- 
bors, and spontaneously fragment. We specified the rates and 
other details about these processes by comparing with many 
different observational constraints. Second, we used the pho- 
tospheric flux sources to extrapolate field lines up into the 
corona by assuming a potential field. 

Despite the model's reliance on flux emergence from be- 
low the solar surface, we did not model the subphotospheric 
motions explicitly. A complete treatment of this problem 
should describe how the photospheric fields are ultimately 



controlled by the overturning dynamics of convection cells 
and their interact ions with one another (e.g-. lFang et alJlToiOl : 
IStein et al]|2010h . In many ways, however, the photosphere 
is believed to act as a relatively "clean" transition layer be- 
tween the highly fragmented fibril fields of the convection 
zone and the space-filling fields o f the corona ( Amariet alJ 
120051: Ivan Ballegooiien & Mackavll20 07j). We take advantage 
of the rapid change in plasma conditions between these re- 
gions to utilize the thin photospheric layer as a natural lower 
boundary. Thus, we used observations of individual features 
and their motions to set up statistical rules for how these fea- 
tures evolve in our Monte Carlo models of the photosphere. 
The ultimate test of the validity of these rules is that the re- 
sulting complex and multi-scale photospheric field matches a 
wide range of observations. (Of course, the observations used 
to test the models must be independent of the observations 
that were used to determine the rules; see Section|4]below for 
more details.) 

Many earUer studies of magnetic flux transport in the 
photosphere were focused on the ne t horizontal diffu- 
sion of fields (e.g., Wang eLal._ 1989i ISimon et alj [1991 
van Ballegooiien et al. 199a). A new era was ushered in by 
iSchrijver et al.i (il997 ). who constructed a statistical model 
that also included flux emergence, cancellation, merging, and 
fragmenting. Numerical simulations of these effects were 
also produced by P arnell (2001), Simon etal. (2001), and 
ICrouch et aP (l2007h . Our Monte Carlo models of the pho- 
tospheric magnetic carpet are based on these earlier models, 
but with three main differences: (1) we use mor e up-to-date 
flux emergence rates (Hagenaar et al.l2()08ll20 1(j). which give 
at least an order of magnitude faster "recycling time" for pho- 
tospheric flux; (2) we model both balanced and imbalanced 
regions on the solar surface that are designed to simulate both 
quiet Sun and coronal hole areas; and (3) we do not presume 
the existence of supergranular motions on the surface — but 
the model does produce a network-like organization of the 
field as a natural output (e.g.. lRasH l2003). 

At each time step in the Monte Carlo simulations, we ex- 
trapolate magnetic field lines up into the corona by assuming 
the field is derivable from a scalar potential. Although the ac- 
tual solar field is likely to h ave signi ficant non-pote ntial com- 
ponents (e.g.. S andman et alj|2009l ; lEdmondson et alj|2009h . 
the approximation of a potential field has been found to be 
useful in identifying the regions where magnetic reconnec - 
tion must be taking place (l Longcopell996l ; IClose et al.ll2005h . 
The potential-field method is also many orders of magnitude 
more computationally efficient than solving the full three- 
dimensional MHD conservation equations. (Doing the lat- 
ter for a system with a complex, evolving, magnetic-carpet- 
like lower boundary is still prohibitively expensive in terms of 
computation time.) Our method involves ignoring the "inter- 
nal" details about how magnetic reconnection actually affects 
the coronal plasma and only investigating the mag netic e nergy 
that is lost via reconnection. We use Longcope's (1 19961) min- 
imum current corona model to take account of the reconnec- 
tion energetics. We emphasize that — despite the title of this 
paper — magnetic reconnection is not a primary "driver" unto 
itself and is merely the end product of the flux emergence, 
cancellation, merging, fragmentation, and diffusion that oc- 
curs on the photospheric lower boundary. 

By modeling only the net changes in the magnetic field 
from one time step to the next, we end up ignoring 
some potentially importan t plasma effects. For example, 
iParnell & GalsgaardI (l2004t) showed that reconnection may 
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progress much more slowly in full MHD than one would ex- 
pect from modeling the syst em as an ide a lized s uccession of 
potential-fi eld states. Als o, 'Lynch et all (|2008|) . iPariat et al] 
y009), Ed mondson et al.l (2010a). and others have shown that 
long-lived, field-aligned currents can exist in the corona due to 
the injection of magnetic flux from below, and these energet- 
ically important structures are not accounted for in potential- 
field models. However, we do not model the most topologi- 
cally complicated regions of the corona, such as the footpoints 
of field lines that connect to the cusps of helmet streamers, or 
to the heliospheric current sheet, or to othe r large-scale sepa- 
ratrix and quasi-separatrix la yers (see, e.g., lEdmondson et alJ 
l2009HAntiochos et al.ll201 0*). Our models generally presume 
the existence of a simple unipolar field at a large height, in 
conjunction with the complex and time-varying magnetic car- 
pet field at the bottom. These "open" unipolar fields may in 
fact close back down onto the solar surface on spatial scales 
larger than our modeled patches of the Sun. Whether this oc- 
curs or not depends on the global distribution of magnetic flux 
across the entire solar surface, which is beyond the scope of 
this paper to model. 

There have been many three-dimensional MHD simula- 
tions of the coronal response to underlying p hotospheric mo- 
tions (see also iG udiksen & Nordlund 2005; Peter et al.ir2006l: 
lGalsgaardll200a llsobe etal.ll2008i) . and this paper does not 
attempt to reproduce those results. The spatial and tempo- 
ral complexity of the footpoint motions in most MHD mod- 
els, however, has usually been assumed to be simpler than 
in the full magnetic carpet as modeled here. We also ignore 
the possibility that there could be a significant back-reaction 
from the c orona on the dynamics of the photospheric foot- 
points (see lGrappin et al.l i2008). Others have studied how the 
evolving photospheric field can affect the properties of coro- 
nal Alfven waves (Malara et al. 2007), coronal mass ejections 
dLvnch et al.ll2009i: iYeates et al. 2010), and the large-scale he- 
liospheric magnetic field (Jiang et al. 2010). The goal of this 
paper is much more limited. We aim to take an initial census 
of the rate at which closed flux opens up from the Sun's mag- 
netic carpet, and to estimate how much magnetic energy may 
be released by the attendant reconnection. Thus, this paper is 
envisioned as a kind of "pathfinder" study that carves out the 
order-of-magnitude expectations for what more sophisticated 
MHD simulations are likely to reveal in detail. 

3. PHOTOSPHERIC FIELD EVOLUTION: MODEL 

In our model, the topology and energy balance of the coro- 
nal magnetic field are assumed to be fully determined by the 
lower boundary conditions at the solar photosphere. Here we 
describe how the photospheric field can be simulated by as- 
suming it consists of a collection of evolving flux sources. 
We developed a FORTRAN code caUed BONES to produce 
Monte Carlo simulations of these flux sources and to trace 
magnetic flux tubes up into the corona. The title BONES was 
inspired by the popular conception of the solar magnetic field 
as a topolog ical skeleton for lo cating important sites of en- 
ergy release dParnell et al.ll2008l) . and also by the dependence 
on randomness in the Monte Carlo technique (i.e., "rolling the 
bones"). 

For a Monte Carlo simulation like this, it is not possible 
to write down a single set of equations that governs the be- 
havior of the magnetic field. Each simulation is a partic- 
ular reaUzation of an ensemble of possible states (see also 
ISchrii ver et al.l [l997l) . Therefore, we must describe the in- 
dividual processes that govern the motion and evolution of 



the flux elements. Section 13.11 introduces some of the gen- 
eral attributes of the BONES simulations. The code models 
the time dependence of the photospheric field as the net result 
of four processes: emergence of new bipoles (Section [3.21 ). 
random horizontal motions (Section 13.31 ), merging and can- 
cellation between pairs of nearby elements (Section lX^l i, and 
spontaneous fragmentation (Section [33] ). 

3.1. Basic Properties and Initial Conditions 

We modeled a patch of the photospheric solar surface as 
a horizontal square box that extends 200 Mm on each side. 
This length scale was chosen to be large enough to encom- 
pass several supergranular network cells, but small enough to 
be applicable to solar wind source regions of roughly uniform 
character (i.e., coronal holes or quiet Sun) and to be able to 
ignore the radial curvature of the solar surface. Thus, the sur- 
face area of the model domain is defined as A = 4 x 10^" cm^, 
or about 0.7% of the Sun's surface area. 

In the part of the BONES code that evolves the photospheric 
magnetic field, each flux element is considered to be a point- 
like monopole having only three attributes: an x position, a y 
position, and a signed magnetic flux $. Even though many el- 
ements are injected into the simulation in equal-and-opposite 
pairs (i.e., as the footpoints of bipole loops), the code retains 
no memory of that association in subsequent time steps. We 
quantized the magnetic flux in units of lO'^ Mx so that incom- 
plete cancellations do not produ ce a huge num ber of infinites- 
imally small elements (see, e.g., Parnell 2001). 

We computed the continuous magnetic field that results 
from the flux elements in several ways. In Section 15.11 we 
describe the computation of the vector field B above the pho- 
tospheric surface. Here we show how an upper limit on the 
magnetic field strength in the flux elements (in the photo- 
sphere) can be used to obtain a lower limit on their spatial 
extent. Let us assume that the horizontal cross section of a 
flux element is circular, and that it is filled with a constant 
vertical magnetic field. It is generally assumed that the field 
in small photospheric concentrations cannot be significantly 
stronger than the so-called equipartition field, in which the 
plasma is in total pressure equilibrium with its (approximately 
field-free) surroundings. In this case, the upper Umit on the 
field s trength is ffmax » 1400 G (see, e.g iParkedl 19761: iLitesI 
I2OO2I: ICranmer & van Ballegooiienll2005h . Thus, we can es- 
timate a lower Umit to the radius of the circular flux element 
as 



l<i>l 



7rB„ 



(1) 



The typical size of observed intergranular G-band bright 
points is re ftj 50-15 km (Muller& Keil 1983). Recently, 
ISanchez Alm eida etal. (2010) measured the filling factor (/ = 
0.89%) and number density (p = 0.97 Mm"^) of bright points 
in quiet Sun regions, and these values are consistent with a 
radius of r^ = (//ttp)'/^ ~ 55 km. The above range of sizes 
corresponds appropriately to fluxes at the low end of the range 
simulated here; i.e., between 10'^ and 10'^ Mx. Elements 
with larger fluxes may not be completely filled by equipar- 
tition fields, and thus they would have larger spatial extents 
than expected from Equation ([T]i- 

At any one time in the simulation, the sum of all positive 
fluxes is denoted $+ and the sum of all negative fluxes is 
denoted $_. These are signed quantities, with $+ > and 
$_ < 0. For all models discussed below that have an im- 
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balance between the two polarities, the sense of the imbal- 
ance is always to have |$+| > |$-|. All results should be 
equivalent for imbalances in the opposite sense. The mean 
magnetic flux densities in the positive and negative flux ele- 
ments, taken over the entire simulation domain, are denoted 
B± = $±/A. Thus, the total "unsigned" or absolute flux den- 
sity is given by Babs = B++\B-\ and the net flux density is 
given by Bnet = |fi+ + fi-| = B+- \B-\. The simulation's flux 
imbalance fraction f is defined as ^ = Bnet/^abs- Small val- 
ues for this ratio (i.e., ^ < 0.3) are typical for quiet Sun 
regions, and larger values (£ > 0.7) are t ypical for coro- 
nal holes ( Wiegelm ann & Sola nki 2004; Z hang et all 120061: 
iHagenaar et aLi2008 : 'Abrame nko et aLi2009D . 

Each run of the BONES code begins with specified initial 
conditions at time f = 0. For models having ^ = 0, there are 
no flux elements in the domain at the beginning of the simula- 
tion. Perfect flux balance is maintained by having all new flux 
elements emerge into the domain at later times as balanced 
bipoles. For models having ^ > 0, the simulation begins with 
a number of identical flux elements, all having positive polar- 
ity, that are distributed randomly over the surface A. These 
initial elements are assumed to each have an equal flux given 
by . 1 times the mean flux in an emerging bipole (see Section 
I3.2l i. The number of these initial elements is determined by 
the input value of the net flux density B„et- As in the ^ = 
case, all new flux elements that enter the domain at f > are 
balanced pairs, and thus Bnet remains exactly constant as a 
function of time. 

For a given simulation that is intended to model a patch of 
the Sun having an imposed flux imbalance ratio ^, the choice 
of the proper input value of Bnet is not known at the outset. 
The overall level of magnetic flux that ends up existing in the 
simulation depends on the collection of dynamical parameters 
that describe the flux emergence, fragmentation, horizontal 
diffusion, and merging (see below). Specifically, the emer- 
gence rate E depends explicitly on ^ (e.g., Hagenaaretd] 
l2008h . Thus, for a given set of dynamical parameters and a 
desired value of ^, we had to produce an iterative set of trial 
runs with a range of guesses for Bnet- Only one unique value 
of Bnet gave rise to a model having the proper self-consistent 
value of ^. After doing this for a range of models, the relation- 
ship between these two parameters was fit with the following 
approximate relation. 



c 



0.268 Bn 



[l + (Bnet/3.58)2-^I] 



2.7110.365 



(2) 



where Bnet > is measured in Gauss and ^ is dimensionless. 
The discrete time step chosen f or the simulatio ns was Af = 
300 s, the same as that used by iParnelll (1200 lb . Five min- 
utes is a representative time scale for photospheric granulation 
(e.g., Ipeubner & Gough 1984), so using a smaller time scale 
would only be appropriate if t he coherent granular rn otions 
were being modeled explicitly. lAsensio RamosI (|2009|) found 
that on spatial scales longer than 300-500 km the solar gran- 
ulation acts as a stochastic, Markovian process. For repre- 
sentative granulation velocities of order 1 km s"' (Hirzberger 
I2OO2I) . this confirms that the minimum resolvable time scale 
(when ignoring coherent convective overturning) should be 
about 300-500 s. For all processes in the BONES code that 
are simulated as occurring st ochastically , we u sed the RAN2 
random number generator of iPress et al.l (Il992h . This routine 
does not repeat its pseudo-random sequence until called at 
least 2 X 10''^ times. This limit was never approached, since 
in even the longest runs of the code the RAN2 routine was 



never called more than lO"* times. 

Over the course of each time step Af , the code updates the 
properties of each of the flux elements from the effects of the 
four general sets of processes described below. 

3.2. Flux Emergence 

Bipolar magnetic features are observed to emerge from 
beneath the photosphere with fluxes spanning several or- 
ders of magnitude from ^10'^ M x (internetwork concen- 
trations) to ^10^^ Mx (sun spots) (ISchrijved 1200 It iParnelll 
2002; Hag enaar et alJ l2008h . Away from active regions, 
much of the emergence tends to occur in the form of bipolar 
ephemeral regions (ERs ) with |$| w lO'^-lO'"^ Mx (see, e.g., 
iHarvey & Martinlll973h . The individual poles of ERs often 
are advected to the edges of supergranular cells and coalesce 
to form network concentrations that end up with similar ab- 
solute fluxes as the ERs themselves (Martin 1988). 

The rate of emergence of ER flux, which we denote E, 
has been estimated in various ways from both measurements 
and models. As the sensitivity and cadence of observations 
has imp roved, the derived emergence rates have generally in- 
creased. ISchriiveJ (1200 lb reviewed earlier measurements and 
models that pointed to a range of E values between about 
2 X 10"^ and 4 x 10~^ Mx cm"^ s"' . Earlier Monte Carlo mod- 
els also found that values in this range seemed to behave in 
similar ways as the real Sun. For example, Parnell (2001J) 
used £■ w 8 X 10"'' Mx cm"^ s~', and Simon et al. (2001) used 
£'wl.3x 10"^ Mx cm"^ s"'. iKrijger & Roudier. (.2003.) found 
that a slightly higher value of 9 x 10"^ Mx cm"^ s"' was 
needed to reproduce TRACE measurements of the chromo- 
spheric network. Assuming a mean flux density in the quiet 
Sun of about 3 to 4 Mx cm"^, it is possible to use the above 
emergence rates to estimate "flux recycling times" between 
about 0.5 and 20 days. 

However, many of these earlier measurements were made 
with sequence s of r elatively low-cadence magnetograms. 
iHagenaar et al.l (^2008) found that when the cadences are re- 
duced from about 90 min to 5 min, many more emergence 
events are obs erved and the emergence rate increases. In fact, 
iMartinI (1198 8*) claimed that it is virtually impossible to even 
identify the same ER from one image to the next unless the 
time cadence betw een them is shor t er tha n about 10 min. The 
revised analysis of IHagenaar et al.l (l2008l) showed that values 
as large as £ w 10"^ Mx cm"^ s~' are often seen in regions 
of balanced magnetic polarities,' along with a noticeable de- 
crease in £■ as ^ increases from to 1. For most values of 
the imbalance ratio (^ < 0.8), these rates of emergence are 
consistent with flux recycling times of only 1-2 hr 

We fit the modified rates shown in Table 2 and Figure 5 of 
IHagenaar et alj ( l2008l 1201 Q) with a quadratic function of the 
imbalance ratio £,, and found 

£ = 7.928 X 10"^(1.356-C^) Mxcm'^s"'. (3) 
For a region with balanced magnetic flux (^ = 0), the max- 
imum value of the emergence rate is £ = 1.075 x lO""* Mx 
cm~^ s~' . As ^ — >^ 1, the parameterized rate declines to a min- 
imum value of £■ = 2.824 x lO"** Mx cm"^ s"'. Note, how- 
ever, that the largest imbalance fraction in the measurements 

' Figure 5 of lHagenaar et alj 120081) showed value s that were erroneously 
reduced in magnitude. The values given in Table 2 of Hagenaar et al.' f200g) 
represented the correct magnitudes for the emergence rates, and a corrected 
revision of their Figure 5 was presented bv Hagenaar et alj >2010iV Our fits 
to these observations utilized a multiplicative correction factor of 5 to the 
numbers shown in their original Figure 5(b), w hich is consistent with the 
updated version shown bv lHagenaar et all 120101) . 
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of iHagenaar et alj ( l2008l) was ^ w 0.94. Our use of values 
larger than this represents extrapolation. It is possible that E 
may decrease more rapidly — possibly to zero — as ^ increases 
from 0.94 to 1 . In any case, we never model the completely 
unipolar case of ^ = 1 . The largest value of ^ used in the mod- 
els presented below is 0.99. 

In order to determine the number of bipoles (A'em) that 
emerge in each time step in the simulation domain, we 
adopted a fiducial value for the average flux per bipole, (<&) = 
9 X 10'** Mx (see below). Thus, A^em =£'AAf/($). In general, 
this does not yield an integer number of bipoles. For a given 
non-integer value of A'em that falls between the two integers n 
and n+ 1, we used the fractional remainder of A'em (in excess 
of n) to determine the statistical chance that the resulting num- 
ber of bipoles is either noin+l. For example, if A'em = 10.22, 
there is a 22% chance that there will be 11 bipoles, and a 78% 
chance there will be 10 bipoles. A new random number is 
generated in each time step to determine whether there will 
be n or n+ 1 new bipoles. 

For each of the emerging bipoles, the BONES code deter- 
mines its total absolute flux by drawing from an empirically 
constrained probability distribution of the form 

P ((b\- i (*-*mi„)eXp [-($-$mi„)/$o] /$i $ > *mm 

(4) 
where the mean flux is given by ( p) = ^min + 2^n . ^ The mea- 
surements shown in Figure 3 of'Hage naaret al.l (l2008 f) pro- 
vided constraints on the functional form of Equation (|4]l, as 
well as values for $min = 2 x 10'** Mx and ($) = 9 x 10"* Mx. 
These values uniquely specify the value of the exponential 
slope $0 = 3.5 x 10'^ Mx. 

In order for the code to sample from the above distribution, 
we computed the cumulative probability distribution by inte- 
grating Equation (|4]l numerically. A parameterized functional 
fit to the inverse of the cumulative distribution was then found 
which allows a uniform random variable (between and 1) to 
be mapped into a proper sampling of Pe(*I')- Once a random 
value of $ has been chosen in this way from the distribution, 
we divided the absolute flux equally between the two poles. 
We note that because the sampling from the distribution is 
random, and because A'em has been truncated to be an integer, 
the exact same amount of flux does not emerge in each time 
step. However, over many time steps the specified emergence 
rate E is maintained on average. 

For each emerging bipole, the x and y positions of the pos- 
itive pole are determined randomly. The position of the neg- 
ative pole is displaced from the positive pole by a horizontal 
distance D and a random orientation angle. The separation 
D must be large enough that the poles will not immediately 
cancel one another out. We assume that D scales with the size 
of the flux element r^ such that D = l.Srcp, where p is the 
dimensionless proximity factor that sets the scale for merging 
and cancellation (see Section [34] |. Since D > rcp, the poles 
are constrained to be noninteracting. For this calculation we 
use the total flux in the entire bipole in the definition of r^- 
(Equation ([Til), so for the mean ($), the mean separation D is 
6.8 Mm. This value of D is within the rather wide observa- 
tional range of separations for newly emerg ed ER bipoles (ap - 
proximately 2-10 Mm), as summarized by iHagenaarl (1200 Ih . 
Note thatliJagenaar (2001) found that D ex $'''^, which is a 
weaker dependence than what we assumed (D ex $°^) by us- 

^ The shape of this distribution is illustrated in Figure|4]below. 



ing Equation ([T]i. 

3.3. Horizontal Motions of Flux Elements 

Magnetic flux concentrations are observed to move around 
on the solar surface in response to plasma flows that occur 
on scales ranging from narrow intergranular lanes (0.05-0. 1 
Mm) up to the supergranular network (~'30 Mm). Our models 
were designed to test the assumption that much of the struc- 
turing on the largest scal es is a natural by-pr oduct of smaller- 
scale motions (see also ICrouch et al.ll2007h . Thus, the mo- 
tions of flux elements are assumed to be of a diffusive charac- 
ter and dominated by granule-scale (1-2 Mm) horizontal step 
sizes. This stands in contrast to other Monte Carlo models of 
the magnetic carpet (e.g., Parnell 2001; Simon et al. 2001) in 
which the motions of the elements are influenced by an im- 
posed supergranular flow pattern. 

For each time step Af, we describe the horizontal motion 
of a flux element as a linear trajectory with speed v and a ran- 
dom orientation angle in the x-y plane. The orientation angle 
is recomputed in each time step with no memory of its pre- 
vious value, so that the long-term trajectory of an element is 
essentially a "random walk." Observationally, the horizontal 
speeds are known to depend on the absolute fluxes in the ele- 
ments, with higher-flux concentrations tending to move with 
lower speeds. Thus, we used a standard exponential fit for the 
mean speed vo. 



VO = Vweak exp 



l$l 



3 X lO'^^Mx 



(5) 



where the constant of 3 x l O'^ Mx in the denomi nator is con- 
sistent with observations (Hagena ar et alj [l999l) and earlier 
models (Sch rijver 200 1). The constant Vweak is the mean speed 
in the limiting case of |<I>| — > 0, and it is a key free parameter 
in these models. The BONES code computes the instanta- 
neous speed V for each flux element by sampling a random 
number from a normal distribution having a mean value of vn 
and a standard deviation of 0.3vo about the mean (see Parnell 
1200 lb . When the horizontal motion is imposed on the x and 
y positions of each flux element, the code assumes periodic 
boundary conditions along the edges of the (200 Mm)^ pho- 
tospheric box. This is designed to take account of elements 
that enter and leave the box via diffusive motions. 

If the horizontal motions were classically diffusive in char- 
acter, the spatial step size Ar could be expressed as 



Ar = V4PAf , 



(6) 



where the diffusion coefficient V is a constant that should 
not depend on the time step Af (see ISchriiveiJl200"Tl) . The 
instantaneous velocity over a single time step would just be 
V = Ar/Af . Solar observations have given rise to a large 
range of values for T), from 50-100 km^ s~' on granular 
scales to 200-2000 km^ s~' on larger scales (e. g.,lBerg er et al.l 
Il998t IHagenaar et al.lfT999l: iGiacalone & Jok ipii 2004). For 
our adopted time step of Af = 300 s, the above range gives 
values of v between about 0.8 and 5 km s"' . 

On granular scales, there is evidence t hat the horizon t al mo- 
tions do not obey classical diffusion. ICadavid et al.l (1 19991) 
found that, for displacement times Af between about 0.1 and 
22 min, the mean-squared displacement Ar^ does not scale 
linearly with Af , but instead 



A^ 



57500 



Af \ 
Iminy 



0.76 



km^ 



(7) 
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For Af = 5 min, this corresponds to an effective velocity 
V « 1.5 km s"'. However, as one examines smaller displace- 
ment times, the instantaneous velocity is larger For Af = 0.1 
min, V increases up to 16.7 km s"'. The observed "subdif- 
fusive" character of the horizontal motions is believed to be 
related to the constraint that flux elements must follow the 
narrow intergranular lanes. Thus, it is not completely valid 
to model the motions as a random walk in a two-dimensional 
plane that ignores the existence of coherent granules. In re- 
ality the elemen ts are constrained to a fractal dimension be- 
tween 1 and 2 (iCadavid et al.l ll999). Even the choice of a 
single value for v may not fully reflect the end-product of un- 
resolved motions taking place within a time step. 

In any case, it is useful to choose a representative value 
for the parameter Vweak that can best reproduce the net disper- 
sal of granule-scale magnetic flux over many time steps. The 
above analysis gives a broad range of plausible choices for 
Vweak between about 0.5 and 20 km s"'. Several trial runs of 
the BONES code were produced with velocities in this range, 
and a final optimized value of v'„eak = 6 km s"' was found 
to produce the most realistic solar conditions. Section |4] dis- 
cusses the results of models constructed with this parameter 
choice. 

3.4. Merging and Cancellation 

In each time step of the simulation, the horizontal distance 
between every unique pair of flux elements is computed. If 
the inter-element distance for a pair is less than a prescribed 
critical value, we assume the flux elements coalesce together 
or cancel one another out. In a computational sense, mergings 
(for like polarities) and cancellations (for opposite polarities) 
are treated in the same way. The flux in the single remain- 
ing element is given by the sum of the two signed fluxes in 
the original elements. The position of this remaining element 
is given by the position of the original element that had the 
larger absolute flux. If an exact cancellation takes place be- 
tween elements with equal and opposite fluxes, then both ele- 
ments are assumed to disappear from the simulation. 

In order to compute the critical distance between a given 
pair of elements, each element is assumed to have a "radius of 
influence" given by rcp, where the constant p is a dimension- 
less proximity factor and r^. is defined in Equation ([T]i. The 
critical distance is the sum of the two radii of influence for a 
pair of elements. 

The proximity factor p is another key free parameter of our 
Monte Carlo simulations. Parnell (2001) essentially assumed 
that p w 2.3 based on an empiric al Gaussian pro file of field 
strength across each flux element. .Schriived (1200 ll) estimated 
the critical mean-free path for interactions between average 
flux concentrations (in quiet network) to be about 4.2 Mm. 
In order to compute a radius of influence consistent with this 
mean separation (i.e., rcP = 2.1 Mm), we can assume that the 
two elements each have a mean flux ($) = 9 x 10'*^ Mx and 
then use Equation ^ to solve for p w 4.6. A series of trial 
runs of the BONES code gave rise to an optimal value of 
p=\Q that produced the most realistic solar conditions (i.e., 
absolute flux densities and number distributions of flux ele- 
ments that agree with the observations discussed in Sections 
4.1-4.2). Thus, for the mean element with ($) =9x 10'^ Mx, 
its radius of influence in the models is 4.5 Mm. 

The BONES code imposes lower and upper limits on the 
radii of influence for the weakest and strongest flux elements, 
respectively. For elements with very low fluxes, the radius of 
influence is not allowed to become smaller than a typical gran- 



ule size of 1 Mm. We assume that the smallest intergranular 
flux tubes can easily traverse the intergranular lan es and inter- 
act in ways that are not resolved explicitly here (iKubo et all 
I2OIOI) . For the strongest flux elements, the radius of influ- 
ence is not allowed to become larger than 10 Mm. Observa- 
tionally, there do not appear to be any mergings or cancel- 
lations that occu r on spatial scales larger than this (see, e.g., 
iLivi et al]|1985h . Practically, though, the imposition of this 
upper limit prevents the occurrence of "long-range" interac- 
tions that would be inconsistent with the existence of the su- 
pergranular network. 

Note that the actual rate of cancellation cannot b e speci- 
fied ex plicitly in these simulations. As described by ' Parnelll 
(I2OOII) . the overall cancellation rate is the eventual result of 
how rapidly the flux elements emerge, move around, and in- 
teract with one another. In a steady state, the cancellation rate 
eventually comes into dynamical equilibrium with the rate of 
emergence E. Thus, our use of the larger values of E from 
iHagenaar et al.l ([2008) implies much mor e rapid cancell ation 
than was found in earher models such as iParnelll (120011) and 
ISimonetal.l(l2001h . 

3.5. Spontaneous Fragmentation 

Observations have shown that magnetic flux elements often 
sphtu p spontaneously into several pieces (e.g.. lBerger & Titi3 
Il996l) . Convective overturning motions on granular scales 
may exert stress on the (usually intergranular) flux elements 
and pull them apart. The physical processes responsible 
for fragmentation are not yet understood, but magnetic re- 
connection may be o ccurring at some stage of the process 
dRvutova et al.ll2003h . There appears to be an observed re- 
lationship betw een the rate of f ragmentation and the total flux 
in an element dSchriiver et al.lll997 ). However, this applies 
only for relatively small concentrations with absolute fluxes 
below about 10^ Mx. Larger concentrations that give rise 
to pores and sunspots tend to survive for longer times, which 
sug gests that the fra gmentation rate saturates for |$| ^ 10^° 
Mx (ISchriiv er"2001'). In our models, we estimated the proba- 
bility of fragmentation Pp (per unit time) to be 



PY'(^)dt 



ko\^dt 



v/i+(|$Mh|)2 



(8) 



where the threshold flux for saturation is given by "I>th = 
3 X 10'^ Mx. This is a slightly simpler version of the param- 
eterization given by Equation (A6) of Schriiver (2001). The 
mean time between fragmentations is given by Pp'. In the 
limit of the largest fluxes, the mean time approaches a con- 
st ant value of (fep^th)"'- 

ISchriiver etaH d 19971) and ISchriiveJ ( 1200 Ih used a combi- 
nation of measurements and models to find values for ko be- 
tween 4 X 10"^^ and 6 x 10"^^ Mx"' s"'. However, these were 
based on the same long-cadence magnetogram observations 
that led to significant underestimates in the emergence rate E 
(see Section [J!2] ). Thus, we decided to increase k^ by approx- 
imately the same relative amount that E was increased from 
the earlier values. The models presented below all use a value 
of/to = 3.5x lO-^-^Mx"' s-i. 

We recompute the probabilities of fragmentation for all flux 
elements in each time step of the BONES code. For cases 
when a uniform-deviate random number (between and 1) is 
less than the probability Pp Af , the code splits the flux element 
into two pieces. The original element keeps a random fraction 
of its original flux (constrained to be between 0.55 and 0.999), 
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and the new element gets the remainder of the flux. The po- 
sition of the original element stays the same, and the new one 
is positioned a distance D away, with a random orientation 
angle. This distance D is the same value discussed in Sec- 
tion l3.2l and it is large enough to prevent subsequent merging 
between the two new flux elements. 

4. PHOTOSPHERIC FIELD EVOLUTION: RESULTS 

In this section we present results from a series of models for 
the photospheric magnetic field as computed by the BONES 
code. A series of tests was first performed to make sure the 
code was actually evolving the flux elements as desired. Once 
the tests verified that each individual process was being mod- 
eled correctly, runs were performed that included all of the 
processes together. We created a basic set of 1 1 models with 
the main adjustable parameter being the flux imbalance ratio 
^. The input values of Bnet for each of these models were iter- 
ated until the final models had steady-state values of ^ equal 
to the desired input values of 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 
0.7, 0.8, 0.9, and 0.99 (see Equation ©). Each model used 
a different integer as a unique seed for the random number 
generator 

As described above, our final Monte Carlo models con- 
tained a much larger emerge nce rate E than did the earlier 
simulations of lParnelll (|2001|) and ISimon eFal] (1200 ll) . If all 
other adjustable parameters had been kept the same as in those 
models, a much larger time-steady magnetic flux would have 
accumulated in the simulation box over time; i.e., the aver- 
aged flux densities would have been much larger than the typ- 
ical values of 3-10 Mx cm"^ observed in quiet regions and 
coronal holes. In order to keep the flux density low, magnetic 
concentrations need to be destroyed as rapidly as they are in- 
jected from below. This is why the BONES code was run 
with more rapid horizontal diffusion (v'weak = 6 km s~'), more 
sensitive merging and cancellation (p = 10), and more rapid 
fragmentation (ko = 3.5 x 10"^^ Mx"' s"') than were used in 
the earlier models. Time will tell if these parameters accu- 
rately represent the real Sun, but as long as the emergence 
rate is high, the models need to facilitate a similarly high rate 
of cancellation in order to produce a realistic steady state. 

Below we present results concerning the overall time- 
stea dy p hotospheric magnetic fields in the simulations (Sec- 
tion |4Tli, the statistical number distributions of flux elements 
(Section 14.2b . and the natural production of supergranular 
magneti c str uctures from the smaller-scale granular motions 
(Sectionl43]l. 

4. 1. General Properties of the Models 

The BONES models were evolved in time, using a step size 
of Af = 300 s, for a total simulation time usually exceeding 
100 days and sometimes exceeding 1000 days (i.e., 10'*-10^ 
time steps). Over the first 10-20 days of a simulation, suffi- 
cient magnetic flux is injected so that the initial conditions are 
completely "forgotten" and the magnetic field reaches a state 
of time-steady dynamic equilibrium. Thus, whenever we cal- 
culate quantities that are meant to represent the time-steady 
parts of a simulation (e.g., means and standard deviations), 
we take only f > 30 days. In the simulated area A, the to- 
tal number of flux elements in the time-steady state tends to 
average between 100 and 200. Although the mean absolute 
flux per injected flux element was ($) /2 = 4.5 x 10'** Mx, the 
eventual mean flux per element in the steady state ended up 
being about a factor of two larger (see below). 



200 




200 



Fig. 1 . — Simulated photospheric magnetograms for random time steps in 
a quiet Sun simulation with § = (top) and a coronal hole simulation with 
^ = 0.8 (bottom). Positive polarities are shown as white, negative polarities 
are shown as black (each saturated at |B;| = 100 G), and the locations of 
magnetic neutral lines (where |B; | = 0) are overplotted as white dotted curves. 

Figure [T] shows simulated magnetogram images for repre- 
sentative time snapshots in two of the models: one for a region 
of balanced magnetic flux (^ = 0) and one for a large degree of 
imbalance (^ = 0.8). The continuous magnetic field strength 
at the photosphere (z = 0) was calculated using the multiple 
monopole model described in Section 15.11 A medium gray 
shade denotes B^ w 0, and the saturation to white and black 
is imposed at B, = +100 and -100 G, respectively. For the 
balanced case, the neutral line meanders through the domain 
stochastically and splits the region into two roughly equal ar- 
eas. For the imbalanced case, the neutral fines surround and 
confine the regions of minority polarity. 

The balanced "quiet Sun" model shown in Figure [Ha) 
has an average total number of flux elements A^ = 163, with 
roughly equal numbers of positive and negative elements and 
an average absolute flux per element of 8.9 x 10'^ Mx. The 
imbalanced "coronal hole" model shown in FigurefUb) has an 
average total A^ = 122, with approximately 8 1 of the elements 
being positive and 41 being negative. Note that if the abso- 
lute flux per element was equal for the positive and negative 
populations, we would have expected that A^( 1 + ^)/2 = 1 1 el- 
ements would be positive, andA^(l-^)/2= 12 elements would 
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Fig. 2. — Time evolution of statistical quantities in the (a) ^ = and (b) 
5 = 0.8 photospheric models. The temporal variability of the box-averaged 
absolute flux density Babs. the total number N of flux elements in the simu- 
lation (divided by 100 to keep the curve in the same plotting domain as the 
other curves), and the flux imbalance ratio § are shown as labeled. 

be negative. Since the number of positive [negative] elements 
is smaller [larger] than predicted, it is clear that the two pop- 
ulations must have different average absolute fluxes. In fact, 
for the ^ = 0.8 model, the average fluxes per element in the 
positive and negative sets were 2.3 x 10'^ Mx and 5.1 x 10'^ 
Mx, respectively. 

In Figure |2] we plot the time dependence of several statisti- 
cal quantities for the ^ = and ^ = 0.8 cases. These models 
reached dynamical equilibrium in only about 5 days of sim- 
ulation time, and only the first 40 days are shown. "* After a 
stochastic steady state has been established, the level of con- 
tinuing tempora l variabihty a ppea rs similar in characte r to the 
simulations of Parnelll (1200 Ih and lCrouchet alJ (l2007h . Note 
that the imbalance ratio ^ does not approach a rigidly constant 
value, but instead fluctuates with a standard deviation that is 
typically 2%-10% of its mean value. 

Comparing Figures|2la) and|2b), we see that as ^ increases 
the mean of the absolute flux density (Bab.s) increases and its 

^ By "dynamical equilibrium" we mean that there appears to be a time- 
steady mean state existing together with substantial variations about that 
mean. It also seems clear that no single ingredient in the photospheiic flux 
evolution model is responsible for determining these time-steady mean prop- 
erties. This state is a complex, nonlinear dynamic balance between emer- 
gence, merging, cancellation, diffusion, and fragmentation. 



variance decreases. Larger values of S, correspond to lower 
rates of flux emergence (see Equation (O), so that a typical 
flux element in the large-^ simulation tends to have a longer 
lifetime before it is destroyed. However, the functional form 
of E{^) is not the only reason for the increase in Babs with in- 
creasing ^. It is possible to illustrate such an increase with 
a simple analytic model that assumes a constant emergence 
rate. If the emergence rate E is fixed, but the box-averaged 
rate of cancellation is assumed to be proportional to the prod- 
uct of the positive and negative flux densities present in the 
box, then their time evolution can be approximated to be a 
simple balance between these two effects, with 

dB+ _ d\B- 

~dt' ~ ~dr 

In a steady state, the time derivatives can be ignored and we 
can solve for E = CB+\B-\. The individual values of the con- 
stants E and C do not need to be specified explicitly, but let us 
assume their ratio E /C is a known constant called Bq. Thus, 
it becomes possible to solve for the absolute flux density in 
closed form. 



E-CBJB- 



(9) 



^abs = B++\B-\ = 



2Bq 



(10) 



The above expression shows how B^bs must increase with an 
increasing imbalance ratio (_, even in the case where E is in- 
dependent of ^. 

Figure [3] shows how the time-steady values of (Babs) from 
the simulations vary as a function of (,. The error bars on 
these model points show ±3 standard deviations around the 
mean values. To ensure that specific realizations of the ran- 
dom number sequences did not affect the results, the means 
and standard deviations for each value of ^ were computed 
from three independent runs of the BONES code. Each run 
used a different random seed, and each run was performed for 
a total of 400 days of simulation time. The modeled abso- 
lute flux densities generally fall between the observationally 
expected limiting values of about 3 and 10 Mx cm"^. Figure 
[3] also shows two curves that illustrate the functional depen- 
dence of the simple analytic estimate of Equation ( fTOb above. 
The two curves, which were computed using the arbitrary nor- 
malization constants Bq =1.4 and 2. 1 G, appear to bracket the 
modeled points surprisingly well. 

In Figure[3]we also plotted measurements made by the Vec- 
tor SpectroMagnetograph (VSM) instrument of the Synoptic 
Optical Long-term Investigations of the Sun (SOLIS) facility 
dKeller et al.ll2003l) . We used publicly available full-disk lon- 
gitudinal magnetograms taken in the Fe I 6301.5 A line. Over 
the time period from August 2003 to November 2009, we ob- 
tained one magnetogram per month for a total of 73 individual 
full-disk maps. For each magnetogram we generated a grid of 
"macropixels" covering the central part of the solar disk (out 
to 0.7 Rq from disk-center). Each macropixel was defined to 
be 100 X 100 raagne togram pixels, or 113" square (see also 
'Hagen aar et alj |2008'). For each macropixel, we measured the 
average flux densities of the positive and negative polarities, 
B+ and B_, and computed Bab,s and ^ as defined in Section lTTl 
A total of 8264 individual measured data points are shown in 
Figure [3] 

The bulk of the low field-strength SOLIS data shown in Fig- 
ure [3] appear to follow the same general increasing trend with 
^ as do the modeled points and analytic curves. The "long 
tail" in the data points that extends upward to 10-100 Mx 
cm"^ represents times when the macropixels covered parts of 
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Fig. 3. — Steady-state dependence of absolute flux density Babs on the flux 
imbalance ratio ^. Mean results from the BONES simulations (filled circles, 
with ±3(7 error bars) are compared with observational data computed from 
SOLIS full-disk magnetograms (gray points), and with the analytic approxi- 
mation given by Equation (To) (dashed curves). 

active regions. Points on the upper-left of the plot represent 
active regions that were mostly centered in the macropixel, 
and points on the upper-right represent times when only one 
dominant polarity of an active region was in the macropixel. 
The models presented in this paper are generally meant to be 
simulations of quiet Sun and coronal hole regions, which are 
sampled by the majority of weak-field data points in the lower 
part of Figure |3] 

4.2. Number Distributions of Flux Elements 

An additional way to verify that the BONES simulations 
produce magnetic fields similar to those on the real Sun is to 
examine the probability distributions of element fluxes and 
compare them with observed distributions. Because the sim- 
ulations typically have only 100-200 elements in them at any 
one time, we sampled the distributions a number of times in 
order to accumulate statistics appropriate for a large number 
of uncorrected patches of Sun. In the models, the time ca- 
dence for this sampling was fixed at 30 days. This time ca- 
dence was found to be more than adequate for the require- 
ment that any given distribution of flux elements must be com- 
pletely recycled from (i.e., uncorrected with) the distribution 
at the previous sampling time. For each case discussed below, 
the simulations were run until the total number of collected 
flux elements exceeded 10^. 

Figure |4] shows example distributions for the two models 
discussed above (f = and 0.8). The distributions of posi- 
tive and negative polarity elements are plotted separately. For 
comparison, the analytic distribution of emerging flux ele- 
ments given by Equation (JUi is also shown. This latter dis- 
tribution has been scaled down in flux by a factor of two 
(i.e., shifted to the left in the plot) to show the distribution 
of fluxes in the individual poles of the emerging bipoles, not 
the total absolute flux in the bipoles as specified by Equation 
(|4|l. For ease of comparison with observations, these p lots are 
shown in the same general for mat as Figures 4 and 6 ofiParnelll 
( 120021) and Figures 2 and 3 of lHagenaar et al.l ( l2008l) . 
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Fig. 4. — Statistical number distributions of flux elements as a function of 
their absolute fluxes in the (a) § = and (b) ^ = 0.8 models. The time-steady 
distributions in the numerical simulations are shown separately for positive 
(solid curves) and negative (dashed curves) polarities, and both are compared 
with the imposed distribution of emerging flux elements (dotted curves). For 
plotting convenience, both the fluxes themselves and the normalized proba- 
bility distributions were divided by lO"*. 

The time-steady distributions shown in Figure |4] are sub- 
stantially "flatter" than the initial distribution of emerging 
flux elements. In other words, the fluxes have spread out 
from the relatively narrow range of injected fluxes (roughly 
lO'^- lO'^ Mx) to both lower and higher values (see Parnell 
I2OO2I) . Most noticeably, the populations of flux elements with 
|<1>| > 3 X 10'^ Mx are hugely enhanced with respect to the 
distribution of injected flux elements. These stronger flux el- 
ements must be the result of mergings between smaller ele- 
ments of like polarity. In addition, the existence of this en- 
hanced strong-flux tail is the reason that the mean flux per 
element is larger than the mean flux in a newly emerged flux 
element (see Section l^TI ). 

Although it is difficult to see in the plots, there is also a 
significant number of elements in the simulations with fluxes 
below the minimum emergent flux per element ($min/2 =10'** 
Mx). These weakest flux elements must be the result of frag- 
mentation and partial cancellation. For the ^ = case, 22% 
of the flux elements have fluxes less than this threshold value. 
Because of their small fluxes, however, these account for only 
about 2.7% of the total absolute flux in the simulation. For 
the f = 0.8 case, 18% of the flux elements have fluxes below 
the emerging threshold value, and they account for 1.2% of 
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the total absolute flux. 

Figure HJb) shows the difference between the distributions 
of positive and negative elements for the imbalanced case of 
^ = 0.8. Overall, the majority polarity has a flatter distribu- 
tion than does the minority polarity, but there is an excess of 
minority polarity elements for the weakest fluxes (|<i>| < 10'"^ 
Mx). This is in good agreeme nt with the observational con- 
clusions of IZhang et al.l (l2006l) for coronal holes. Also, the 
differences in shape shown in FigureUb) are highly reminis- 
cent of the flux element distributions shown in Figure 2 of 
iHagenaar et alj (i2008 ) for coronal holes. 

4.3. Naturally Occurring Supergranular Scales 

The resemblance between the cellular pattern of solar 
granulation and that of the larger-scale supergranulation 
has long been interpreted as evidence that both phenom- 
ena are manifestations of the Sun's convective instabil- 
ity (e. g., 'Leighton et al." '1962"; 'Roxburgh & TavakoH 119791: 
ISimon & Weiss 1991; Rieutord & Rincon 2010). However, 
because the flow patterns in the supergranular network are 
weak and intermittent, it has not been possible to definitively 
prove their convective origin. It may be that multiple interac- 
tions between granule-scale structures produce a distributed 
network of downflows that in turn seeds horizontal super- 
granular flows and the aggre gation of strong network fields 
(lRastJl2003l: iGoldbaum et alJl2009.) . Alternately, the opposite 
may be the case; i.e., it may be the aggregation of small-scale 
magnetic fields that gives rise to the weak supergranular flows 
dCrouch et al. 2007). In this section we show that the BONES 
simulations provide some evidence for the initial magnetic- 
field aggregation described in the latter scenario. 

How are the spatial scales of supergranulation measured? 
It is well known that the dominant cell sizes are of order 
10-30 Mm, but different types of measurement give differ- 
ent answers. Simon & LeightonI ( 11964 ) found cell diame- 
ters around 32 Mm by interpreting a utocorrelation functions 
of chromospheric Dopplergrams. ISingh & Bappul (Il981l) 
traced the cells manually, based on Ca 11 K-line intensity im- 
ages, and f o und diameters of ^22 Mm. Wang (1988) and 



IWang et aP (1199 6) applied the autocorrelation technique to 
magnetograms and found scale sizes between 10 and 25 Mm, 
depending on the precise d iagn ostic techniques used. Finally, 
iDe Rosa & Toomrd (12004 !) and'Ha genaaret all (Il997h used a 
range of sophisticated algorithms to trace and characterize su- 
pergranular boundaries, and found average diameters of only 
-15 Mm. 

Because the BONES simulations predict only the proper- 
ties of the magnetic field — and neither the chromospheric 
emission nor the Doppler velocities — we decided that the 
most straightforward comparison to make would be with the 
measu red magnetogram autocorrelation functions of Wan g 
(Il988l) . First, a random time step from each of the 1 1 mod- 
els was used to create simulated magnetograms similar to 
those shown in Figure [T] Then, for each y row in the mag- 
netogram, we computed a series of one-dimensional autocor- 
relation functions in the x direction for the scalar value of B,, 
i.e., 

/+00 
B,(x, y) B,ix+x',y)dx, (11) 

oo 

which was then normalized such that AC(0,3') = 1. Figure|5ja) 
shows an example autocorrelation function from the ^ = sim- 
ulation, plotted as a function of the lag parameter x' . Similar 
results were found when the roles of the x and y coordinates 
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Fig. 5. — (a) Example of a simulated magnetogram autocorrelation func- 
tion for a slice across the g = model, plotted as a function of the lag param- 
eter x' (see Equation jilt ), (b) Results for modeled mean values of FWHM 
(filled circles) and SM (open circles) plotted as a function of ^, with error 
bars denoting ila in the simulated dist ribut ions of values, and the observed 
ranges of FWHM and SM values from iWang U98S) (gray regions). 

were reversed. 

We characterized the model autocorrelation functions by 
finding both the full-width at half-maximum (FWHM) of the 
central peak and the distance between the central peak and the 
next secondary maximum (SM). Doing this for each value of 
y gave rise to ensembles of values for FWHM and SM in each 
of the 11 simulations. Figure|3b) shows the mean values for 
each of these ensembles, along with error bars that show ± 1 
standard deviations about the means. There is no significant 
^ dependence in the modeled values. For all 1 1 simulations, 
the average model FWHM is 4.48 Mm and the average SM 
distance is 25.1 Mm. These va lues compare favorably to the 
solar observations reported by 'Wanj] (Il988h (shown as gray 
bars in Figure |5]), who found FWHM values between 4 and 6 
Mm, and SM distances of 15 to 20 Mm. 

The benefit of making a direct comparison between sim- 
ulated and observed FWHM and SM values is that there is 
no need to interpret these quantities in terms of arbitrarily 
defined cell diameters."* The models appear to succeed in 
roughly reproducing the observed autocorrelation properties 
of the network. It may be possible to explain this success by 

■* See, however, Figure[8]below for a more intuitive way of visualizing the 
naturally occurring "supergranular network" in these simulations. 
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invoking processes of diffusio n-limited aggregation as sug- 
gested by ICrouch et alJ (l2007h . In this picture, time-steady 
magnetic structures "collect" on specific scales that depend 
on the combined emergence, diffusion, and cancellation of 
flux elements. Supergranular flows may the n occur as a result 
of the magnetic structuring. Crouch et al.l (12007) performed 
tests with a Monte Carlo model that varied several of the dis- 
crete step sizes and interaction distances, and found that the 
resulting supergranular scale size does not depend on these 
input parameter choices. Instead, it is the overall level of flux 
emergence and horizontal diffusion — which in turn drives the 
cancellation rate — that sets the time-steady distance between 
network concentrations. 

5. CORONAL FIELD EVOLUTION 

One of the major goals of this paper is to explore how the 
complex photospheric fields in the magnetic carpet connect 
with time-variable open flux tubes and closed loops in the ex- 
tended corona. Thus, here we describe how the field lines are 
traced upwards and are evolved in time (Section lSTt . we sum- 
marize the resulting open and closed fields as a function of the 
flux imbalance ratio ^ (Section l572] i. we compute relevant time 
scales for the opening up of closed flux tubes (Section I5.3l l, 
we estimate the amount of magn etic energy that emerges in 
the form of bipoles (Section 15.4b . and we compare it to the 
energy released into the solar wind by magnetic reconnection 
(Section [ 



5.1. Field-line Extrapolation Method 

As summarized in Section |2] we compute the vector mag- 
netic field B above the photospheric surface by assuming the 
field is derivable from a scalar potential. In other words, each 
flux element is assumed to act as a monopole-type source, 
with 

B« = E?/^' (12) 

where the coordinates r, = (xi,yi,Zi) specify the locations of 
each flux element /, and the field point r = ix,y,z) can be lo- 
cated anywhere at or above the photosphere (z > 0). ^ is the 
flux in each element (see, e.g., IWang|IT998t 



signed magnetic fl 
ICloseetalJ '20Q3). 



To avoid singularities at the solar surface, all ele ments are 
assumed to be "submerged" below the photosphere dSeehafeiJ 
119861: iLongcopa i2005i) . For simplicity we assumed that all 
flux elements are at a constant depth. We chose an optimum 
value of Zi = -1 Mm on the basis of the following consid- 
erations. The peak magnetic field strength Bpeak in the pho- 
tosphere, due to a single flux element, occurs right over the 
point itself at x = x,, y = y,, and z = 0. Thus, 



B 



peak ~ 



2nzj 



(13) 



We want to ensure that |Bpeak| is less than the equipartition 
field strength Bmax for all elements in the simulation (see Sec- 
tion ITT)- Because we do not model pores and sunspots, we 
can apply this constraint to elements up to a maximum flux of 
I'J'I w 10^*' Mx. Thus, applying the condition |Bpeak| < Bmax to 
Equation (fT3] l for this value of the flux gives rise to |z,| > 1.1 
Mm. On the other hand, observations have shown that the 
field strengt h in a recentl y emerged ER is at least a few hun- 
dred Gauss (Martin 1988). For the average flux in one pole of 
an emerging ER (i.e., ($)/2 w 4.5 x 10"^ Mx), we apply the 
condition Bpeak ^ 100 G and obtain an upper limit |z,| < 0.85 



Mm. The two above constraints on the magnitude of z, are 
formally incompatible with one another, but the value ^ 1 Mm 
appears to be a likely compromise between the two. 

The BONES code contains a subroutine that can either trace 
field lines up from the photospheric surface or down from a 
larger height. The incremental path length As for numerical 
steps taken along the field varies with height, from a minimum 
value of 0.03 Mm at the photosphere to a maximum value of 
10 Mm at a height of z = 200 Mm. At intermediate heights. 

As = (0.03 Mm)'-'^(10 Mm)'^ , (14) 

where C = z/(200 Mm). Field lines that begin at the photo- 
sphere are traced until they either curve back down to intersect 
the z = plane again (and are called "closed") or they climb 
past a maximum height of 200 Mm (and are called "open"). 
As discussed in Section |2l on the real Sun it is possible that 
many flux tubes that reach higher than 200 Mm may eventu- 
ally be closed back down in the form of large-scale helmet 
streamers. Whether this occurs or not depends on the global 
distribution of magnetic flux across the entire solar surface. 
In any case, it is likely that some plasma that reaches large 
height s in streamers also interacts with the accelerating solar 
wind ( Wang et al.l 1200 0*). so it may not be too erroneous to 
classify these field lines as open. 

When the Monte Carlo simulation of the photospheric field 
settles into a dynamical steady state (defined here as f > 50 
days), we begin tracing field lines in order to compute the 
coronal vector field in each time step. This essentially as- 
sumes that any temporal changes occur "instantaneously;" 
i.e., with a time scale shorter than Af = 5 min. In similar kinds 
of potential-field simulations, Regnier (2009) found that the 
actual delay between a given photospheric impulse and the 
response higher up in the corona is only of order 2 min. Thus, 
our assumption that B(r) can be recomputed from each time 
step's new lower boundary condition appears to be reason- 
able. 

In order to quantify the changes that occur in the magnetic 
field from one time step to the next, we trace a set of field lines 
that is associated with the A^ flux elements on the surface. The 
general idea is to compare the open/closed topology of flux 
tubes that can be identified unambiguousl y both at the begin - 
ning of a time step and at the end (see also lClose et al.ll2()05l) . 
If a flux element moves around on the surface and does not 
undergo substantial merging, cancellation, or fragmentation, 
then we can say that it has "survived" that time step, and thus 
it makes sense to evaluate how its open/closed connectivity 
may have changed. In cases where the merging, cancellation, 
or fragmentation makes only a minor change to an original 
element's flux, we also consider that element to have survived 
when the element's flux changes by less than a specified frac- 
tional threshold S. In most runs of the BONES code presented 
below, (5 = 0.1. This means that if a flux element ends the time 
step with a flux that is within 10% of its original flux, it is clas- 
sified as being the same element. Flux elements that cannot 
be tagged in this way are not counted. We discuss the effects 
of varying the S parameter below. 

Rather than just trace one field line from each flux element, 
we instead chose to more finely resolve the coronal magnetic 
field by tracing seven field lines from each element. The 
initial footpoints of these seven field lines are arranged in a 
hexagonal pattern with respect to each flux element's circu- 
lar "patch" on the surface. One field line is centered on the 
flux element. The other six are arranged in a ring around the 
central point with an angular separation of 60°, each at a hori- 
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zontal distance of rc{l+p)/2 from the central point. This dis- 
tance is halfway between the flux element's intrinsic radi us re 
and its critical interaction distance as defined in Section [3741 
At the beginning of each time step the BONES code traces 
IN field lines and tags each footpoint with a unique (nonzero) 
numerical identifier. Each of the flux tubes associated with 
element i is assigned an equal magnetic flux $,/7. During 
the progress of each time step, new flux elements that emerge 
are given an identifier of zero. Also, if merging, cancellation, 
or fragmentation changes the flux in an element to a degree 
greater than the relative threshold S, its numerical identifier is 
reset to zero. At the end of each time step, the coronal field is 
traced again for the subset of surviving flux elements that have 
nonzero numerical identifiers. The magnetic flux in those el- 
ements is grouped into four bins that are defined by whether 
the flux tubes were open or closed at the beginning of the time 
step, and whether they are open or closed at the end. Section 
I5.2l discusses the distributions of magnetic flux in those four 
bins. 

We note that our method of accounting for the open and 
closed magnetic flux has several potential shortcomings. By 
not counting either the newly emerged flux elements or those 
that undergo substantial merging, cancellation, or fragmenta- 
tion, we run the risk of not seeing fields that may be releasing 
lots of energy via magnetic reconnection. We will see below, 
though, that the magnetic-carpet evolution is not so vigorous 
that these flux elements represent a significant fraction of the 
total number. In fact, for most models the fraction of magnetic 
flux that is missed by not counting these "rapid evolvers" is 
only of order 5% to 15%. Another possible limitation of our 
method is that we trace the identities of individual flux tubes 
for only one time step. If we wanted to measure more ac- 
curate time scales for flux reconfiguration, it may have been 
advantageous to follow field lines for more than just one time 
step. However, since the magnetic carpet keeps evolving, the 
number of flux tubes that would become uncountable (i.e., 
missed by virtue of exceeding the threshold 5) increases for 
each additional time step over which flux-tube survival would 
be traced. Following field lines only over the course of one 
time step, with Af = 5 min, gave the best balance of time res- 
olution and flux capturing. 

5.2. General Results 

Figure |6] illustrates a selection of field lines for BONES 
models with a mostly balanced lower boundary (^ = 0.2) 
and a highly imbalanced lower boundary (£_ = 0.8). The 
three-dimensional field lines are shown projected into a two- 
dimensional plane that is defined by an observer viewing the 
scene at an inclination angle 82° from the normal to the photo- 
sphere. Two different shades denote closed versus open field 
lines. Models with more imbalanced fields (i.e., higher values 
of have both a larger fraction of open flux and a smaller 
vertical extent for the closed loops. Both of these trends are 
examined quantitatively below. 

We studied the statistical properties of the closed loops in 
the simulations by tracing large numbers of field lines from 
random starting locations {x,y,0) in the photosphere. Exam- 
ple time snapshots from the 1 1 models (with varying ^ values) 
were used to trace at least 5000 loops in each model. For the 
six models with ^ < 0.5, for which there were fewer open 
field lines, we were able to compute at least 20000 loops. The 
maximum heights of these loops were collected into 1 1 statis- 
tical distributions, one for each model. Although the means 
and standard deviations of these distributions were computed. 




(a)e = 0.2 




(b)? 



Fig. 6. — Traced magnetic field lines at example time steps in BONES 
models having (a) ^ = 0.2, and (b) ^ = 0.8. Open and closed field lines are 
plotted in black and gray, respectively. In both panels, the horizontal box 
outlines the (200 Mm)^ photospheric simulation domain. The vertical scaling 
has been stretched by about a factor of two, such that the uppermost tips of 
the field lines are at a height of z~ 110 Mm. 

the distributions were far from Gaussian in shape. Thus, we 
quantified them further by computing percentile intervals H„ 
of the sorted cumulative distributions of heights. For example, 
25% of the loops have heights less than the quartile height of 
H25, and 50% of the loops have heights less than the median 
height of HsQ. We also computed Hjs and Hg^, with the latter 
being an approximate indicator of the largest loops (without 
being dependent on the statistically insignificant tail of the 
very largest loops). 

Figure I2] shows how the percentile intervals vary as a func- 
tion of the flux imbalance ratio ^. On the smallest spatial 
scales (i.e., for granule-sized loops characterized by H2S and 
7/50) there does not appear to be a significant dependence on 
^. However, the longest loops follow the trend that is visually 
apparent in Figure|6j i.e., the more balanced the photospheric 
field, the larger the loops. This trend is apparent not only in 
7/75 and 7/95, but also in the mean height (//) that is weighted 
more strongly by the longest loops. 

Figure |7] also shows approximate observational ranges 
of mean loop heights for quiet Sun (QS) and coronal 
hole (CH) regions as determined by 'W iegelmann & Solankil 
(2004). These loop-height calculations were similar to 
ours in that they were based on potential-field extrapo- 
lations from photospheric low er boundary conditions, but 
IWiegelmann & Solankil ( |2004|) used observed magnetograms 
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Fig. 7. — Variation of percentile intervals of the sorted statistical distri- 
butions of loop heights, shown as a function of ^. Percentiles at the 25%, 
50%, 75%, and 95% levels (solid curves) are compared with the mean loop 
height (H) (dashed cu rve) and with observationally inferred values from 
IWieeelmann & Solankil <2004') for quiet Sun (QS) and coronal hole (CH) re- 
gions (gray boxes). 

from the Miche lson Doppler Imager (MDI) instmment on 
SOHO (see also IClose et ani2003t man et alJl2010l: llto et alJ 
I2OO 8). The overall agreement with the modeled ^ depen- 
dence of (H) is good. The general trend for high-^ CH re- 
gions to have shorter loops than low- g QS regions is also 
consi stent with the trend p ointed out by lFeldman et alJ (Il999h 
and lGloeckler et al] (I2003D for the source regions of fast solar 
wind to be correlated with short loops and the source regions 
of slow wind to be correlated with long loops. 

A representative illustration of the footpoints of open field 
lines is given in Figure [8] for the ^ = 0.8 model. This plot 
shows the locations of the photospheric footpoints of IC* field 
lines that were traced down from an evenly spaced grid at 
the top (z = 200 Mm). In order to account for the horizon- 
tal flaring of potential field lines from the finite-sized simula- 
tion box, the grid of 100 x 100 starting points had an overall 
horizontal size of 1800 x 1800 Mm in the x and y directions 
(centered on the 200 x 200 Mm simulation box). The overall 
appearance of Figure [8] is highly reminiscent of the observed 
supergranular network. The apparent "cell diameters" tend to 
be between 20 and 40 Mm as on the real Sun. Note also the 
appearance of thin channels, stretched between smaller knots 
of closed-field regions, that appear to support the connectivity 
theorems described by lAntiochos et al.l (l2007h . 

All of the 10"* open field lines with footpoints shown in Fig- 
ure [8] are of positive polarity. This is the dominant polarity as 
specifie d by the initial conditions of the BONES code (see 
Section lTTT i. All negative polarities end up connected to pos- 
itive polarities in closed loops, and thus there are no "open 
funnels" with the non-dominant polarity. Of course, this is 
also a highly simplified situation when compared to the real 
Sun, for which there are often network concentrations of both 
polarities even in strongly unipolar coronal holes. 

As described above, at the beginning of each time step there 
is a set of field lines traced from each of the flux elements. 
These IN field lines are used to estimate the instantaneous 
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Fig. 8. — Photospheric locations of footpoints of "open" magnetic field 
lines traced down from an evenly spaced grid at a height of j = 200 Mm, for 
one time snapshot of the ^ = 0.8 model. 

fractions of absolute unsigned flux that are either open or 
closed. The fraction of flux that is open is denoted /open, and 
in Figure|9]we show its mean value as a function of the £_ im- 
balance ratio. This fraction is never exactly the same from 
one time step to the next, and the error bars show ± 1 stan- 
dard deviations about the mean values. On average, /open is 
roughly equal to ^ itself. In other words, models with bal- 
anced fields tend not to have much open flux, but when there 
is an increase in the unbalanced component of the field there is 
a corresponding increase in the fraction of open flux. Figure|9] 
also compares the modeled values of /open with observational 
determinations of this quantity from IWiegeknann & Solankil 
(2004), and there is a similar trend of direct proportionahty, 

with /open « £.■ 



5.3. Comparison of Relevant Time Scales 

We studied the time evolution of magnetic topology in the 
BONES simulations by following the opening and closing of 
flux tubes from the beginning to the end of each time step. 
For comparison, we also computed the recycling time scale 
for flux to emerge from below the photospheric surface (see 
also Section [372] l. We defined this quantity as 



{Babs) 



(15) 



For our models we took (Babs) from Figure |3] and E from 
Equation ([3]), and we found that the emerg ence time scale 
Tem te nds to have values around 1-2 hr (see iHagenaar et alj 
I2OO8I) . Regions of extreme flux imbalance undergo slower 
emergence, with Tem exceeding 10-20 hr when ^ > 0.9. Fig- 
ure fTOf a) shows the ^ dependence of this time scale. 

Next we used the flux tubes traced in our simulations to in- 
vestigate the_timescalesfor magnetic field evolution in the 
corona. IClose et aP (l2005h performed a similar study in the 
limit of a balanced field, with ^ = 0. They computed a so- 
called coronal flux recycling time that is meant to character- 
ize a local rate of change of the coronal field. This rate is 
driven both by reconnection and by topological evolution of 
the complex "hierarchical tree" of footpoint domains in the 
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Fig. 9. — Various dimensionless flux fractions shown as a function of t^: 
mean values of /open (filled circles) and mean values of »/) (open circles), both 
with their JiIct spreads shown as error bars, and observational estimates of 
/open from lWiegelmann & SolankH 120041) in QS (triangles) and CH (squares) 



magnetic carpet. Because changes in tiie coronal field can 
take place even w ithout any flux emergence or cancellation, 
IClose et al] (120051) found that coronal flux recycling times 
can be significantly shorter than photospheric flux recycling 
times. Changes in topological connections can occur purely 
as a result of the horizontal m ot ions of flux elements (e.g., 
lEdmondson eTani2009l l2010ah . IClose et al.l (l2005h used an 
older photospheric flux recycling time of T^m ~ 15 hr, but they 
found that the coronal flux recycling time can be as short as 
1.4 hr. When emergence and cancellation were suppressed, 
the coronal time scale was approximately a factor of two 
larger (~3 hr) but still m uch more rapid th an Tsm- Our mod- 
els differ from those of Close et al.l (l2005h in that our pho- 
tospheric emergence time scale is now of the same order of 
magnitude as their coronal recycling time scale. 

Below we describe how we estimate how long it takes for 
just the open flux to recycle itself in the corona. We do not 
track the (possibly more numerous) changes in topology that 
do not involve open flux tubes. As summarized in Section lSTTl 
over the course of a time step some of the flux in the model is 
unaccounted for because it has either emerged since the last 
time step or it has evolved beyond recognition as the same flux 
element. The remaining fraction of total absolute flux — i.e., 
that which survives the time step unaltered — is called ifj, and 
Figure |9] shows how its mean value increases steadily from 
about 0.82 to 0.95 as ^ increases from to 1. A larger choice 
for the relative tolerance parameter 5 would give a larger sur- 
vival fraction if) (see below), but it can be argued that too much 
tolerance would give rise to errors in how flux tubes are iden- 
tified and tracked. 

For flux tubes that survive a time step relatively unchanged, 
we compared the endpoints of the field lines traced at the be- 
ginning and end of the time step. The fluxes in these field lines 
are summed into four separate bins that are defined by their 
connectivity. The four bins correspond to four fractions of the 
total surviving absolute flux: foo (starts open, ends open), /oc 
(starts open, ends closed), /co (starts closed, ends open), and 



/cc (starts closed, ends closed). Because the overall magnetic 
configuration of the system does not vary strongly over a sin- 
gle time step, we found that foo ~ /open- Also, the two frac- 
tions that denote change (/oc and /^o) both tend to be small 
contributors to the total. The mean values of /co in the models 
tend to vary between about 0.005 and 0.025, with the largest 
values occurring for intermediate imbalance ratios of ^ ss 0.5 
and the smallest values occurring at the extremes of <^ = and 
0.99. We also note that the time averages of /co and /oc are 
always roughly equal to one another (as should be required 
for a time-steady dynamical equilibrium). For all 1 1 models, 
the time averages of these two fractions never differ from one 
another by more than about 2%. 

At any one time, we define the amount of open (absolute) 
flux density as Bopen = /openBabs- We computed the instanta- 
neous rate of opening in each time step Af as 



dB 

dt 



/co^al 
Af 



(16) 



Note that the above definition makes the implicit assumption 
that /co is the fraction of the total absolute flux density in the 
simulation that opens up in one time step. However, this frac- 
tion is only approximately i/; times the total absolute flux that 
opens up. We assumed that the small fraction (1 - V') that was 
not counted contributes in the same way as the larger fraction 
if) that was counted. (This assumption is tested below.) Thus, 
the mean time scale for the opening up of closed flux tubes is 
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(/co) 



(17) 



Because the quantities /co and {dB/dt)^^ can be quite vari- 
able from time step to time step, we realized that care should 
be taken in computing the averages in Equation ( fTTb . We 
ended up computing these averages in two independent ways. 
First, we took simple arithmetic averages of the time series 
for (dB/dt)co and the other quantities. Second, we integrated 
the rate defined in Equation ( fT6l ) as a function of time to build 
up the cumulative amount of flux density that is opened up 
over the course of the simulation. This is a monotonically 
increasing function, but its increase with time is intermittent 
because different amounts of flux are opened up in each time 
step. We fit the cumulative growth of opened flux density with 
a linear function, and then used the slope of this linear fit as 
the mean value of {dB/dt)^o- These two methods gave results 
that agreed with one another to within about 10%, and we 
used the latter technique for all values reported below. 

Figure [TOl' a) compares the above time scales with one an- 
other. It is clear that Tco ~ r^m in these models. In other words, 
the time scale for the replacement of the photospheric flux — 
via emergence from below — is the same as the time scale for 
replacement of the open flux that feeds the solar wind. At first 
glance, this appears to be a simple requirement for a time- 
steady equilibrium, in the same way that /co ~ /oc is required 
to maintain a steady state. However, one can imagine situ- 
ations where the rate of flux evolution in the corona is not 
so strongly coupled to the em ergence rate of new flux from 
below (e.g.. lClose et alJl2005h . In our case, it is the use of po- 
tential fields — which are remapped during each time step with 
no allowance for the storage of free energy in the corona — 
that demands Tco ~ r^m- In other words, the BONES models 
reproduce the case of highly efficient magnetic reconnection, 
where the corona "processes" the flux as quickly as it is driven 
(stressed or injected) from below. One can imagine that in 
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" Standai'd value used in all other models discussed 


below. 











Fig. 10. — Comparison of time scales for various models, (a) For the Monte Carlo models of the magnetic carpet, the recycling time for flux emergence (dotted 
curve) is compared with the time scale for flux opening (filled circles and solid curve), (b) For the Cranmer et al. 1 2007) solar wind models, we plot acceleration 
times T„i„i] up to heights of 25 Mm (dashed curve), 50 Mm (dotted curve), 100 Mm (solid curve), and 200 Mm (dot-dashed curve) versus the outflow speeds at 1 
AU. Also shown is an approximate region of parameter space that corresponds to upper heights z, that exceed 2-3 times the maximum heights of closed loops in 
the corresponding BONES models (gray box). 

a full MHD simulation the efficiency of magnetic reconnec- 
tion may not be so high, and thus the resulting non-potential, 
current-filled corona should exhibit t^o > Tem- 

Note that Figure [TOf a) does not show the value of t^o for 
the ^ = model. As Equation ([TtI i makes clear, in this case 
both the numerator and denominator are numbers that should 
approach zero. Ideally, there should be no open fields at all in 
a perfectly balanced potential field. The BONES models do 
in fact give slightly nonzero values for (/open) and (/co), but 
these are believed to be numerical artifacts arising from the 
discrete nature of the field-line tracing technique. We reiterate 
that we do not compute the time scale for all of the coronal 
flux to be recycled. That recyclin g time should be n onzero 
even for the balanced ^ = model (IClose et al.ll2005h . In all 
models with ^ ^ 0, the full coronal recycling time is likely to 
be significantly shorter than t^o- 

In order to study the dependence of our results on the as- 
sumptions made about flux-tube identification, we varied the 
threshold flux identification parameter 6 away from its default 
value of 0.1, in a range between and 0.5. This parame- 
ter sets the relative tolerance for the classification of evolving 
flux elements over a time step. Table 1 shows several result- 
ing parameters of the test simulations, which were all per- 
formed for ^ = 0.4. As we expected, the flux survival fraction 
ip increases monotonically with increasing 6. However, there 
does not seem to be any definitive trend with 5 in the frac- 
tion of flux that opens up (/co), the related time scale for flux 
opening (tco), or the energy flux relea sed b y reconnection into 
open-field regions ((Fco), see Section l53T l. This suggests that 
the topological changes resulting from flux-tube opening are 
adequately resolved in the simulations. Thus, we retain the 
standard value 6 = 0.1 for the remainder of the paper 

It is worthwhile to compare the time scale for flux opening 
to the time scale for solar wind acceleration along the open 
flux tubes. If a significant amount of solar wind plasma flows 
out during the time it takes the open field to reorganize itself 
via reconnection, then the reconnection processes themselves 
probably are not responsible for producing the majority of 
the solar wind. The RLO idea depends on the plasma in open 
flux tubes coming from the opening up of closed loops. Thus, 



we want to determine whether or not a large amount of mass 
accelerates out in the open flux tubes over the time it would 
take for significant mass to be processed via loop-opening. 

The time scale for wind acceleration from a lower height 
Ztr in the solar transition region (TR) to an arbitrary upper 
height z is 

Twi„d(z) = / 4^ , (18) 

^ZTR «(Z ) 

where uiz) is the radial wind speed. The TR was chosen as 
the height to start the integration because that is where the 
mass flux of the wind is thought to be determined (see , e.g., 
iHammedfigsllWithbrodfTosllHan steen & Leei^'1995'). We 
used the one-fluid solar wind models of Cranmer et al. (200'$ 
to compute Tv^ind, and we defined zxr as the height at which 
the temperature in a wind model first reaches 10^ K. 

Figure fTOlb) shows the wind acceleration time scales for 
several representative upper heights z, and for a range of mod- 
els of the fast and slow wind that have speeds at 1 A U be- 
tween 400 and 750 km s~' (see ICranmer et al. 2007). Two 
side-by-side plots are necessary in Figure [lO] because there 
is not a unique one-to-one correspondence between the flux 
imbalance ratio ^ and the wind speed at 1 AU. We do know, 
however, that there is some association between slow wind 
streams and QS regions on the surface (^ w 0) and between 
fast wind streams and CH regions on the surface (^ « 1). 
Thus, the overall left-to-right variations in the two panels can 
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be roughly associated with one another. 

The slow solar wind models shown in Figure [TOl b) have 
the shortest acceleration time scales. Given Equation ( fTSl l. 
this is potentially counter intuitive. However, we note that the 
slow wind models from Cranmer et al.l (12007 ) often have lo- 
cal maxima in u{z) of order 100 km s~' in the low corona 
that are not present in the mor e steadily accele rating fast wind 
models (see also Figure 7a of Cranmenl2010l) . These regions 
correspond to enhanced magnetic fields that were included to 
simulate open fields at the edges of streamers and active re- 
gions. Observati ons are beginning to show hints of such rapid 
outflows as well (iHarra et al.ll2008l: ISubramanian et al.ll2010t 
iBryans et aLlbOlOh ! 

When comparing the time scales for flux opening and solar 
wind acceleration, we can use the loop heights illustrated in 
Figure |7] as an order-of-magnitude guide for the maximum 
height z to use when computing Twind(z)- For example, when 
parcels of solar wind exceed a height that is 2-3 times Hgs, it 
can be safely assumed that the wind has left behind virtually 
all interactions with closed loops and should be considered 
to be freely accelerating. This allows us to compare the time 
scales between panels in Figure [TO] for the two general types 
of solar wind: 

1 . For slow wind streams rooted in balanced QS regions 
(i.e., ^ w 0), the height at which the wind flows "free 
and clear" of loops is of order 50-100 Mm. Figure 
[TOl b) shows that this height corresponds to Twind ~ 0.1- 
0.3 hr This is a shorter time scale than the representa- 
tive flux-opening time Tco ~ 1 hr that corresponds to the 
left side of Figure [TOf a), but it is still of the same order 
of magnitude. Thus, it is possible that RLO processes 
could be important for slow wind acceleration. 

2. For fast wind streams rooted in unbalanced CH regions 
(i.e., ^ « 1), the height corresponding to 2-3 times //gg 
is only of order 5-15 Mm. The fast wind accelerates to 
this range of heights in less than about 0.3 hr, but the 
flux-opening recycling time in coronal holes can be as 
long as 3-10 hr. This is a larger discrepancy than in 
the case of the slow wind, and it implies that it is un- 
likely that RLO processes are important in accelerating 
the bulk of the fast wind. (Of course, it still may be 
the case that RLO processes produce a highly intermit- 
tent or episodic injection of mass and energy into the 
fast wind in coronal holes — just not enough to affect 
the majority of the accelerating plasma. The polar jets 
discussed further in Section|6]may be a prime example 
of this intermittency.) 

The gray box in Figure [TOl b) shows the approximate range 
of wind acceleration time scales that correspond to maximum 
heights z exceeding about 2-3 times Hys as discussed above. 
The shape of the gray region is roughly independent of wind 
speed and ^. This is because, as one goes from left to right 
in the plot, the increase in Twind (for constant z) is offset by 
the fact that the relevant value of z decreases (because Hgg 
decreases; see Figure |7]l. 

Finally, we reiterate that the values of t^o shown in Figure 
[TOl' a) are likely to just be lower limits to the actual time scales 
of flux-opening. As discussed above, our models assume a 
succession of potential fields that are consistent with the as- 
sumption of rapid magnetic reconnection. If the true MHD 
state of the corona exhibits less efficient magnetic reconnec- 
tion, then the photospheric footpoint stressing will build up 



non-potential fields and current sheets in the corona and thus 
give rise to larger net values of Tco- In this case, it is even 
more certain that t^o ^ Twind, and our conclusion that RLO 
processes are unimportant in accelerating the solar wind is 
strengthened. 

5.4. Poynting Flux in Emerging Bipoles 

Our primary reason for constructing the BONES simula- 
tions was to estimate how much energy is deposited into the 
solar wind by the evolving magnetic carpet. First, though, it 
is necessary to compute how much magnetic energy is being 
injected into the system from below the photosphere. It is 
not obvious that all (or even most) of this energy is able to 
be converted into forms that supply heat or momentum to the 
accelerating solar wind. Since, on small scales, much of the 
injected magnetic energy is in the form of compact bipoles, it 
may be difficult for much of this energy to become "liberated" 
into the open-field regions when these bipoles evolve and in- 
teract with one another Thus, in this section we discuss the 
total magnetic energy that is potentially available, and in the 
following section we estimate what fraction of it is actually 
released by reconnection into the open-field regions. 

The relevant quantity to compute when considering the rate 
of injection of magnetic energy from below the photosphere 
is the Poynting flux, which is defined as 



S=^ExB«--L[(vxB)xB], 
47r 47r 



(19) 



and where the latter approximation assumes the ideal condi- 
tion of MHD flux freezing. In the Cartesian system studied in 
this paper, the most relevant component of the Poynting flux 
is the z component, with 



S,= —[B\v,-{y^-'Ri_)B:] 



4tt 



(20) 



where Bj^ and vj^ are the components of the magnetic field 
and velocity in the horizontal (x-y) plane. The two terms on 
the right-hand side of Equation ( l20t represent components as- 
sociated with flux emergence and surface flows, respectively. 
For simplicity, though, in the remainder of this section we 
will endeavor only to estimate the overall magnitude S of the 
Poynting flux. This gives a reliable upper limit that is inde- 
pendent of the adopted geometry and topology of the emerg- 
ing flux elements. 

Observationally, the Poynting flux can be estimated from 
various measured proxies (e.g., Welsch et al. 2009), but there 
exist ambig uities in the da ta that give rise to significant uncer- 
tainties. .Fisk et al.l (Il999h estimated the magnitude of S to be 
about 5 X 10^ erg cm~^ ^~ ' i" sou rce regions of the solar wind. 
iMartfnez Gonzalez et al.l (1201 Oi) used vector magnetic fields 
measured by Hinode/SOT to estimate that small-scale emerg- 
ing loops provide something like 10^ to 2 x 10^ erg cm"^ s"' 
to the low chromosphere in quiet regions. 

We estimated the magnitude of the Poynting flux for the 
Monte Carlo models developed above in two independent 
ways. Figure [TTT a) shows that the two methods gave rise to 
similar ranges of Poynting flux (both of order 10'' erg cm"^ 
s"') with a relatively weak dependence on (,. These two meth- 
ods are described below. 

First, we note that the emergence rate E (Equation (O) al- 
ready describes how much magnetic flux is driven up from 
below the photosphere, per unit area and per unit time (i.e., its 
units are Mx cm"^ s"'). What we want to know is how much 
magnetic energy emerges, in units of erg cm~^ s"'. Thus, if 
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Fig. 1 1. — Compai'ison of energy fluxes for various models, (a) Estimated flux (Fco) in loop-opening events (filled circles and solid curves) computed with 
two choices for SlCl. Also shown are approximate Poynti ng fl uxes 5 for photospheric flux emergence, with the dotted region showing estimates from Equation 
)23t and the gray region showi ng estimates from Equ ation J26). The dashed curve shows a linear scaling (Fco) c< ^. (b) Total dissipated solar wind energy flux 
F„jnd from the WTD models of lCranmer et alj <2007D . 

is the height-dependence of B, for z > 0, that is the major 
source of uncertainty in eval uating Equ ation ( l24l i. However, 
it is straightforward to follow iFisk et alJ ( 1999) and assume a 
vertical falloff that depends on a power of heliocentric radius. 
Thus, 

B^Bq (^) (25) 



we can relate the flux in an emerging bipole to its magnetic 
energy, we can convert easily from E to S. Treating a pair of 
equal-and-opposite emerging flux elements as an an ideal (but 
partially submerged) magnetic dipole, we can specify its field 
strength as 

:^Vl+3cos20, 



B = 



(21) 



where $,■ is the absolute flux in each pole, D is the horizontal 
separation between the two poles, r is the distance measured 
from the center of the dipole, and 6 is the polar angle mea- 
sured from the (horizontal) dipole axis. Assuming the dipole 
is submerged at a depth |z,|, it is possible to integrate the mag- 
netic energy L'mag over the full coronal volume V (i.e., over all 
X and y, and all z > 0) analytically. We thus found 



(where r = z+Rq), and then 



BqRq 



2n2 



<^fD 



f/mag- I dV ^^- ^28^2U.|3 



(22) 



87rTein(2«-l) 

At large distances above the photosphere, the exponent n ap- 
proaches a value of 2, but it is be lieved to take on larger va l- 
ues closer to the surface (see, e.g.. lBanaszkiewicz et al.ll998h . 
For a typical value of Bq = 4 G and rem = 1 hr, we can esti- 
mate an upper limit on S by assuming n = 2, and thus obtain 



Note that the magnetic energy above the photosphere is ex- 
tremely sensitive to the submerged depth |z,|. Once the mag- 
netic energy due to a given bipole is known, we can estimate 
the magnitude of the Poynting flux as 



For a more realistic coronal value of 

l5 „,„ 0,^-2 c-1 



{u. 



($) 



(23) 



where the angle brackets denote the properties of the "av- 
erage" emerging bipole as discussed in Section 13.21 Figure 
[TIT a) shows this quantity for the 1 1 models as a function of ^, 
and for two reasonable choices of |z,| (0.8 and 1.2 Mm). For 
typical values of £ = 10"^ Mx cm"^ s"^ ((())= 9 x 10'^ Mx, 
D = 6.8 Mm, and \zi\ = 1 Mm, we find that 5 « 8 x 10^ erg 

-2 -I 

cm s 



The second way to estimate S was proposed by iFisk et all 
( Il999l) . Here, we compute the total magnetic energy in the 
system (per unit surface area) and divide it by the flux recy- 
cling time. In other words. 



1 f B^ 
S^— dz—. 

Tern . on 



(24) 



Here, the value of B at the photospheric surface is essentially 
the time-averaged absolute flux density (i.e., Bq sa (Babs))- It 



S = 4-x lO^ergcm'^s"'. 

n w 8, we have 5 w 6 x 10^ erg cm"^ s"'. Figure [TTT a) shows 
how S varies as a function of ^ when the modeled variations 
in (Babs) and rem are used, and when the two above values of 
n = 2 and 8 are assumed to define the lower and upper limit- 
ing cases. Given the uncertainties, the two alternate methods 
of estimating S give numerical values that are quite consistent 
with one another 

5.5. Energy Release in Loop-Opening Events 

We used the output of the BONES simulations to estimate 
the amount of energy released by magnetic reconnection for 
cases of closed flux tubes turning into open flux tubes (and 
vice versa). It is important to note that there are also expected 
to be many other sites of reconnection and energy release 
that do not involve open flux tubes. For example, in a bal- 
anced QS region there may be a large number of small-scale 
"footpoint-swapping" events that start with a configuration of 
closed loops and end wit h a slightly different topological dis - 
ti-ibution of closed loops (iPriest et al.l2002HClose et al.l2005l) . 
In this paper, we explicitly ignore the energy release in the 
closed-closed events in order to focus on only the subset of 
events that can input mass and energy into the solar wind. 

The basic geometrical picture for a flux-opening event is 
the "anemone" type structure that is believed to exist at the 
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footpoints of many X-ray bright points, coronal jets, and po- 
lar plumes (e.g. , Svrovatskii 1982; Shib ataet al.lll992ll2007l: 
iFilippov et ar[l2009t iShimoio & Tsunetall2009h In this pic- 
ture, a small bipolar magnetic field either emerges or ad- 
vects into the presence of a larger-scale open field. Mag- 
netic reconnection is believed to occur roughly above the 
end of the bipole with the oppo site polarity as the open field 
jEdmondson et al.ll2009ll2010al) . The newly opened flux may 
take the form of a jet or plume (jWang 1998), and the newly 
closed flux may "subduct" and provide heating to the underly- 
ing chromosphere (Gugli elmino et al. 2008). In one of these 
interchange-reconnection type events, the amount of closed 
magnetic flux that opens up should be the same as the amount 
of pre-existing open flux that becomes closed (i.e., /co ~ foe)- 
Because we model the evolution of the coronal magnetic 
field as a succession of potential fields (see Section 13, we 
use the quasi-static "minimum current corona" (M CC) model 
to estimate the energy loss du e to rec onnection ('Longcop^ 



1996HLongcope & Kankelborg 1999; Be veridge & Longcope 



2006"). In this model, the mean energy flux released in 
closed-to-open reconnection events is proportional to the rate 
{dB/dt)co at which magnetic flux is opened up (see Equation 
(fT6b). For our simulations, we derived the MCC energy flux 
to be 



9lCl 



{d) 



dB 

dt 



(27) 



where (E>i is the mean absolute flux per element, {d) is the 
mean separation between elements in the simulation, and 6'l 
and Cl are dimensionless constants. The Appendix presents 
a detailed derivation of Equation ( |27| i for anemone-type re- 
connection events, including a discussion of the most likely 
numerical values for 6'l and Cl. 

It is important to clarify that the energy flux given by Equa- 
tion dZTJ l is meant to be an order-of-magnitude representation 
of the magnetic "free energy" released by reconnection. The 
MCC model depends on an estimate of the current that builds 
up and is dissipated along an idealized separator, and truly 
non-potential MHD simulations are needed to verify whether 
these estimates are valid. Also, the MCC model does not 
specify how the energy is partitioned into other forms such 
as thermal energy, bulk kinetic energy, waves, MHD turbu- 
lence, and energetic particles. Determining this partitioning 
is a complex problem — one definitely beyond the scope of 
this paper — that often requires the use of fully kinetic simu- 
lations. However, it has been found that many forms of parti- 
cle energization that occur rapidly and locally in reconnection 
regions may eventually become unstab le to dissipation that 
randomizes the velocity distributions (iBhattacharied 120041; 
iFuiimoto & Machidal2006t lYamadal2007 '). Thus, much of the 
energy that initially goes into, e.g., waves or supra-Alfvenic 
beams may end up released in the form of heat. This will 
be our implicit assumption when comparing F^o with the en- 
ergy fluxes required to heat the corona and accelerate the solar 
wind along open flux tubes. 

Figure [iTT a) shows the time-averaged quantities {Fco) for 
10 of the standard BONES models as a function of ^ (ex- 
cluding the case ^ = 0). See below for a discussion of how 
Fco varies in time. The lower and upper sets of points were 
computed by assuming the product of the two dimensionless 
constants ^lCl to be 0.003 and 0.011, respectively (see the 
Appendix). For nearly all of the models, (Fco) is significantly 
smaller than the available Poynting flux S. For the lowest val- 
ues of ^, the resulting "efficiency" of energy release in open- 



field regions (i.e., {F^q) /S) may be as low as 0.001-0.01. This 
means that in QS regions, only a tiny fraction of the magnetic 
energy that enters the system ends up being available for driv- 
ing the solar wind via RLO processes. 

For most values of ^, the computed values of (Fco) are sig- 
nificantly lower than the canonical heat fluxes ( i.e., 3 x 10^ to 
10'' erg cm"^ s"') that lWithbroe & Novesid 19771) estimated are 
needed to maintain QS and CH regions on the Sun. However, 
for the most unbalanced CH regions (^ > 0.95) the modeled 
energy fluxes do appear to approach both the empirically re- 
quired heating rates and the empirically constrained Poynting 
fluxes. Observed coronal holes, ho wever, exhibit values of 
^ over a mu ch wider range of values (IWiegelmann & Solankil 
2004; Abra menko et al.l2009h . so the models still have a prob- 
lem with explaining CH coronal heating in general. 

Figure [TTT a) also shows a curve that represents a linear de- 
pendence with the flux imbalance ratio; i.e., (Fco) oc ^. For 
0.2 < ^ < 0.9, this linear relationship appears to fit the varia- 
tion in the modeled energy fluxes. Because we also know that 
/open oc S, (see Figure |9]l, this tells us that the heating rate in 
flux-opening events is roughly proportional to how much of 
the time-averaged magnetic field remains open. 

As was done in Section |53] above, we can also compare 
the results from the BONES simulations with earlier models 
of solar wind acceleration along open flux tubes. We would 
like to assess how much energy flux needs to be deposited 
in open-field regions in order to produce t he solar wind. 
We us ed the one-fluid WTD-type models of ICranmer et alJ 
(l2007l) to estimate this quantity. These models involved find- 
ing a self-consistent description of the volumetric heating rate 
Q= I V • F| (in units of erg cm"-* s"') that was able to maintain 
time-steady corona and solar wind. In order to derive the to- 
tal energy flux |F| that was dissipated in one of these models, 
we had to integrate over the entire radial grid, which extended 
from the photosphere to the heliosphere. The ICraruner et alJ 
(2007) models were computed along magnetic flux tubes that 
have a radially varying cross-sectional area Atube(z)- Thus, the 
radial integral of the product QAtube gives the total power dis- 
sipated (in erg s~') in a flux tube. To express this quantity as 
an energy flux and compare it to the quantities shown in Fig- 
ure fTTT a). we normalized the area function Atube(z) to the area 
of the simulation box (A = [200 Mm] ^) at a height correspond- 
ing to the low corona, at which the supergranul ar funnels have 
expand ed to fill the "canopy" volume. For the lCranmer et alj 
( |2007|) models, this height corresponds to z sa Q.QARq ~ 28 
Mm. Then the energy flux can be computed by dividing the 
total power by the box area A, and 

^wind = |F| = T / dzA,^Uz)Q{z) . (28) 



FigurefTTTb) shows how F„ind depends on the wind speed at 1 
AU for the same models that were shown in Figure [TOl' b). We 
point out that .Fisk et al.l (11999) was correct to conclude that 
the energy flux needed to accelerate the solar wind is of the 
same order of magnitude as the emerging Poynting flux S (see 
also Leer et al. 1982; Schwadron & McComas 2003). How- 
ever, Figure [m a) shows that RLO-type flux-opening events 
do not appear to be able to release the required energy flux 
into the open flux tubes. 

A key result of many c oronal heating models — including 
the MCC models of lLongc ope ( 1996) — is that the energy dis- 
sipation process should be highly intermittent. This occurs 
in the BONES simulations as well. Figure [T2l shows a snap- 
shot of the time dependence of the quantity Fco for the ^ = 0.2 
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Fig. 12. — Time evolution of the energy flux Fco released by reconnection 
into open-field regions, for BONES models having ^ = 0.2 (gray curves) and 
^ = 0.8 (black curves). Time averages for both cases are denoted by dashed 
lines. 

and 0.8 models. These heating rates were computed with the 
upper-limit value of the product O^Ch = 0.01 1 . Thus, the time 
averages of these quantities correspond to the upper set of 
solid points in Figure [TTT a). For the majority of the mod- 
els (0.2 < ^ < 0.9) the standard deviation of Fco is approxi- 
mately half of its mean value. For the extreme models with 
the lowest and highest values of ^, the standard deviations in- 
crease to be about equal to their means. Such a scaling would 
be expected if the energy fluxes were sampled from an ex- 
ponential distribution similar in form to that of the emerging 
bipole fluxes (Equation (|4|i). In any case, the variability of the 
predicted heating rates may be just as useful as the mean val- 
ues when attempting to d istinguish between different coronal 
heating models (see, e .g., [ Parked 1 9881: f Walsh & Galtieil2000t 
iBuchUn & Venill2007h . 

It is worthwhile to list some of the ways in which the above 
models may be incomplete or incorrect. For example, 

1 . The assumption of a succession of potential fields is 
likely to limit the verisimilitude of the models. It 
is clear that time dependent, three-dimensional MHD 
models — which contain currents, resistivity, and finite- 
pressure effects — would shed more light on the dynam- 
ics and energetics of this system. If the gas pressure 
in localized reconnection regions begins to exceed the 
magnetic pressure (i.e., /3 ^ 1), there may be additional 
ways for the flux tubes to "break open" that were not 
accounted for here. 

2. Even within the confines of a succession of poten- 
tial fields, the assumptions of the MCC model may 
be too simplistic. For example, it is known that in 
three-dimensional reconnection there are both spatial 
and temporal variations of the current along separators, 
whic h our implementation of MCC does not incl ude 
(e.g.. [Galsgaard & Parnellll2005l; IParneU et al.ll2010l) . 

3. Our assumption of ^l = 1 in Equation (|27] | may be too 
large, and thus our resulting estimate for the energy flux 
released by reconnection may be too high. 



4. The simple three-pole magnetic geometry discussed in 
the Appendix did not consider realistic asymmetries 
in either the footpoint locations or the magnitudes of 
the flux sources. When such asymmetries are included 
(lAl-Hach ami & PontinI 1201 Oh . the resulting range of 
values for Cl would Ukely be different. It is unclear 
whether Cl would be larger or smaller than the values 
estimated in the Appendix. 

5. The use of the mean flux element separation (d) in 
Equation (|27] | is only a rough approximation. Since 
there may be significant energy release when one flux 
element gets very close to another, it may be better to 
use a mean distance that is smaller than (d). In that 
case, our estimate for the heating rate could be too low. 

6. As we mentioned in Section 15.11 above, many of the 
flux tubes that are classified as "open" may in fact be 
closed in the form of hydrostatic helmet streamers. In 
reality, then, the energy flux that escapes out into the 
solar wind could be even lower than the values of (Fco) 
that were shown in Figure fTTT a). It is also possible 
that large-scale interchange re connection c ould even 
tually open up these flux tubes dWang et al.1^2 000; Fis 
120051: lAntiochos et al.ll2010HEdmondson et al.ii2010 ' 
but modeUng those processes is beyond the scope of 
this paper. 

Roughly speaking, there appear to be just as many reasons 
why our results for the rates of RLO heating and flux-opening 
may be overestimates as there are reasons why they may be 
underestimates. Despite the approximate nature of these mod- 
els, however, we believe that the main result (i.e., {F^o) <C S 
for most values of ^) is not likely to be wrong by many orders 
of magnitude. 

6. DISCUSSION AND CONCLUSIONS 

The primary aim of this paper was to begin testing the con- 
jecture that the opening up of closed flux in the Sun's mag- 
netic carpet is responsible for driving the solar wind. First, 
we created Monte Carlo simulations of the complex photo- 
spheric sources of the solar magnetic field. The resulting 
time-averaged properties of the models appeared to agree well 
not only with observations of the flux density and the flux im- 
balance ratio, but also with observed probability distributions 
for the flux elements and autocorrelation functions of the field 
strength. A supergranular pattern of network magnetic con- 
centrations appeared spontaneously in the models, despite the 
lack of any imposed supergranular motions. Then, armed with 
some degree of confidence that the model photosphere is an 
adequate reflection of reality, we then computed the coronal 
magnetic field. Assuming that the coronal field evolves as 
a succession of potential-field extrapolations, we were able 
to estimate both the time scales and energy fluxes associated 
with RLO-type flux-opening events. 

From the simulations, we found that the Poynting flux in 
emerging magnetic elements (which could be a proxy for the 
maximum energy flux available for coronal heating) is typ- 
ically around 10* erg cm"^ s"'. However, for quiet regions 
((, ^ 1), only a tiny fraction of the available Poynting flux 
was found to be released in flux-opening events via magnetic 
reconnection. A similar situation was found to exist in mixed- 
polarity regions that can correspond to either quiet Sun or 
coronal holes (^ < 0.8). For the most unbalanced coronal hole 
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regions (^ sa 1), the fraction of Poynting flux released in flux- 
opening events may approach unity. In these regions, how- 
ever, the time scale for flux opening was found to be signifi- 
cantly longer than the solar wind travel time from the coronal 
base to heights far above the tops of loops. Thus, it appears 
that a significant amount of mass accelerates out into the so- 
lar wind over the time that it would take for the plasma to be 
processed via RLO type mechanisms. From the above esti- 
mates of time scales and MCC energetics, we conclude that it 
is unlikely that the solar wind is driven by reconnection and 
loop-opening processes in the magnetic carpet. 

Despite the negative conclusion regarding the solar wind 
as a whole, we believe that the physical processes mod- 
eled in this paper are likely to be relevant in many other 
ways. For example, it is possible that more can be learned 
about the energetics of polar jets with the methodology de- 
veloped here. Soft X-ray observations can be used to esti- 
mate the energy flux released due to jet eruptions. These jets 
are believed to span several orders of magnitude in the to- 
tal amount of energy released ; i.e., betw een about 1 0^^ and 
10^'^ erg (Shimoio et al. 1998; iChifor et al. 2008 ; Pariat et alj 
120091; lMoritaetalJl2010t) . Let us tak e a canonical value o f 
Ejet « 4 X 10f_ergjrom the model of ' Shimojo eTaP (Il998l) . 
Recently. ISavcheva et al.l (2007) identified 104 jets with the 
Hinode X-Ray Telescope (XRT) over a time span of 44 hours 
in a polar coronal hole, which gives a mean time between jets 



(for the observed area ) of • 



jet' 



1500 s. The area examined by 



ISavcheva et al.l (12007 ) was approximately the "front half" of 
the polar cap, viewed from the side, which extended down to 
about 25° colatitude and thus covered about Ajet « 1.5 x 10^' 
cm^. Thus, we estimate the mean energy flux released in jets 
to be i^et ~ £jet/(^jetTiet) ~ 2 X 10"* crg cm"^ s"' . This agrees 
reasonably well with the predicted energy fluxes (for ^ w 0.6- 
0.9) shown in Figure fTTT a). 

The flux-opening events modeled in this paper may also be 
releva nt to under standing the small eruptions seen in quiet re- 
gions (llnnes et a l. 2009; Schrijver 2010) that may be related 
to coronal mass ejections (CMEs). However, it is not guar- 
anteed that every jet-like eruption observed in the corona re- 
leases material that accelerates up into the solar wind. There 
is observational evidence that at least some coronal jets con- 
tain plasma tha t falls back down because it failed to re ach the 
escape speed ( Baker et al.l 120081; IScuUion et al.ll20()9l) . This 
may put some jets into the same category as spicules, which 
are known to carry orders of magnitude more mass up (and 
down) than is needed to feed the solar wind (e.g.. Sterling 
HOOO; De Pontieu et al. 2009). 

A potentially valuable set of observational diagnostics of 
the processes discussed in this paper are the elemental abun- 
dances and ionization states of different particle species that 
escape into the solar wind (Zurbuchen 2007). The closed-to- 
open reconnection events that we have modeled may inject 
some plasma with a distinctly "closed" composition signa- 
ture into regions that have signatures otherwise dominated by 
flux tubes that remain open. It is worth noting, however, that 
there remains disagreement about exactly what kinds of abun- 
dance and ionization sign atures signal the pres ence of closed 
loops, and which do not. Cranmer et al.l (l2007l) showed that a 
range of WTD-type open-flux-tube models can produce val- 
ues of the commonly measured O'^IO^^ and Fe/O ratios that 
agree reasonably well with in situ measurements (see also 
iPucci et al.l 12010 ). Thus, we question the popular assertion 
that the charge-state and first-ionization-potential (FIP) prop- 



erties measured in the slow solar wind can only be explained 
by the injection of plasma from closed-field regions on the 
Sun. 

Whether or not the solar wind energy budget is accounted 
for by RLO processes, the inherent variability in the mag- 
netic carpet is likely to cause some kind of MHD fluctu- 
ations to propagate up into the corona. The response of 
the coronal field to the evolving footpoints may result in 
Alfven waves with peri o ds of order Tem ~ Tco (see Figure 
[T0|. In fact, iHoUwed (Il990l |2008) suggested that "flux 
cancellation" events in the corona may be the most likely 
source of the long-period (i.e., 0.5-10 hour) Alfven waves 
that dominate in situ measurements. The statistical proper- 
ties of these low-frequency fluctuations may also be consis- 
tent with an ori gin in the motions of c oronal fiel d-line foot- 
points ( Matthaeus & GoldsteinI Il986t iGiacalone & Jokipiil 
l2004tlNicoletal.ll2009l) . 

In order to further test the applicability of RLO- 
type processes to accelerating the solar wind, the mod- 
els need to evolve beyond the approximate potential- 
field "skeleton" and to incorporate MHD effects. Multi- 
dimensional MHD simulations (e.g ., JGudiksen & Nordlunc 
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2008t [Edmondson et al. 2009) illustrate the aspects of coro- 
nal reconnection that are — and are not — modeled well by po- 
tential fields, and future studies need to account for these ef- 
fects more consistently. Also, analytic models of the micro- 
scale kinetic physics should be developed further in order to 
complement the coarser-gridded numerical simulations. Id eas 
such as stochastic growth theory (Cairns & Robins onll998h or 
non-modal stability (Camporealeet al. 2010) may be useful 
ways to understand the partitioning of energy within recon- 
nection regions. 

Additional work should be done to refine and test the idea 
that the supergranular network is the natural by-product of 
smaller-scale granular activity (Rast 2003). Our success in 
reproducing the measured autocorrelation patterns in magne- 
tograms (see Figure|5]l does not necessarily imply that there is 
no convective component to supergranulation. However, our 
results do appear to provide evidence that at least some of the 
10-30 Mm magnetic structure on the Sun can be built up from 
^1 Mm granulatio n effects via a kind of diffusion-limited ag- 
gregation (see also lCrouch et alJl2007h . 

Another topic that requires further study is the cou- 
pling between waves and flux emergence in the gran- 
ular convective flows at the ph otospheric lower bound- 
ary. ICranmer & van BallegooijenI ({2005) estimated that the 
surface-averaged energy flux of Alfven waves in the low 
corona is of order 10^ erg cm~ ^ s~' (see Figure 12 of 
ICranmer & van Ballegooijenll2005h . It is probably not a co- 
incidence that this is of the same order of magnitude as the 
Poynting flux S due to the emergence of ephemeral regions. 
The interplay between convective overturning motions, col- 
liding granular cells, and thin flux tubes may give rise to a 
rough equipartition between these different sources of energy. 
By constructing models that contain the seeds of both WTD 
and RLO processes, we can better determine their relative 
contributions to coronal heating and solar wind acceleration. 
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APPENDIX 
AN IDEALIZED APPLICATION OF LONGCOPE'S MCC MODEL FOR ANEMONE-TYPE EVENTS 

In this section we show how the 'Longcogg d 19961) MCC model can be applied to the results of the BONES simulations 
described above. In this model, the motions of discrete flux sources on the solar surface give rise to stresses in the coronal field 
that are concentrated at topological boundaries (i.e., separatrix surfaces and separator field lines). Electric currents are assumed 
to form along the separators, and then dissipate as magnetic reconnection occurs in response to the evolution of the flux domains. 
iLongcopa (Il996h found that the power dissipated in a single flux transfer event must be choppy and intermittent, but its time 
average can be written as 

- „ /* <i$ 



'"^Yc 



dt 



(Al) 



where d^/dt is the time derivative of magnetic flux that is in the process of transferring its connectivity, /* is a characteristic 
current that is assumed to flow along the separator, 0l is a dimensionless threshold constant describing the intermittency of 

reconnection, and c is the speed of light i n vacuum . 

In the double-bipole configuration of Longcope ( 1996), the transfer of magnetic flux (d^/dt) occurred because a fraction of 
the flux from the positive pole of one bipole became reconnected with the negative pole of the other bipole. In our model, we 
consid er the transfer of flux from a closed flux tube to an open flux tube, or from open to closed. Equa tion (B9) of iLongcopg 
(1 19961) gave the characteristic current used in Equation dAlb above. Correcting a typographical error in lLongcopd (Il996h . this 
current is given by 

'■-'-^ 

where L is the length of the separator field line, s is a dimensionless geometrical constant (with s=\ corresponding to a circularly 
shaped separator field Une), and B^ is an average value of the Jacobian-like perpendicular derivative of the vector field at the 
separator, 

Bl = V-det(Vj_Bj_) . (A3) 

In the above, the perpendicular direction is defined relative to the separator field line. 

For a given magnetic configuration, the above equations let us estimate the power emitted from the loss of magnetic free energy 
via reconnection. However, it would be too computationally intensive to locate and trace all of the separator field lines during 
every time step of the BONES simulation. Thus, we aim to simplify the application of Equation (lAlb by creating a characteristic 
"building block" for the magnetic geometry in a typical (anemone-type) opening/closing event. These building blocks can then 
be assembled together in a statistical way to account for the total amount of evolving flux in each time step of the Monte Carlo 
simulations. 

Wang ( 19981) described a simple model of plume/jet events in coronal holes that involved only three discrete flux sources: two 
that form a localized bipole and a third that represents a unipolar source of open field. As discussed in Section 1531 the energy 
release that is assumed to occur in this system happens when some of the flux in the bipole reconnects with the unipolar region, 
giving rise to an equal amount of opening and closing of flux (/co = foe)- For geometric simplicity, let us assume that all three 
flux sources are collinear along the x axis, with a negative source in between two positive sources. The flux evolution occurs as 
the negative pole of the bipole moves away from its original positive partner and towards the positive source of open field. We 
want to evaluate the properties of this system at a representative time in the middle of its evolution, so let us posit an additional 
symmetry; i.e., we assume that the negative pole sits at the origin (x = 0) and the two positive poles are both equidistant from the 
origin (x = -izd) and of equal positive flux. This may be an extreme simplification, since it is known that many details of three- 
dimensional null-point reconnection do depend on whether the geometry is symmetric or asymmetric (lAl-Hachami & PontinI 
I2OIO). However, the other uncertainties in the order-of-magnitude MCC model are probably not outweighed by this issue. 

To evaluate the coronal magnetic field arising from this three-pole system, we set the flux in the positive poles to $+ > and 
flux in the negative pole at the origin to $_ < 0. The two free parameters that constrain the topology of the field lines are the pole 
separation d and the ratio of negative to positive fluxes m = \^-/^+\. Thus, Equation (fT2l i gives 

$+ f x+d x-d mx 1 

.M,y,z) = — I [(_^^ ^^2 + 3,2 + ^2] 3/2 + [(_^_^)2 + 3,2 + 22]3/2 " [^l+yl + ^ip/l j ' ^ ■> 



^+ ( y y my 

2^ \ [(jC + £/)2 + 3;2 + z2]3/2 "^ [(x-£/)2 + / + z2]3/2 " [x2 +y2 + ^2j3/2 



By(x,y,z) = — <j r.,. . .^2 . ,.2 . -7i3/2 + r.,. .^2 . ,.2 . -2iV7 " r,.2 . ,.2 . .0,^/0 I ' (^5) 



^^r z z mz \ 

z{X,y,Z) - — I [(^_^^)2+3;2 + ^2]3/2 + [(_^_^)2 + 3;2 + ^2]3/2 [-^2 +y2 + ^2j3/2 j ' ^ > 

We will consider values of the flux ratio m between about 0.5 and 2. For m > 2, the central source "breaks out" with its own open 
field of negative polarity, which is a situation that we are not considering here. 

Figure [T3l a) illustrates a few representative field lines for the case m = 0.8. For simplicity, we assume the poles are at z = 0. 
The coronal volume (z > 0) can be separated into four distinct domains according to the field-line topology: (1) a set of open field 



22 



CRANMER AND VAN BALLEGOOIJEN 




10 



10" 



10 
1.0 

0.1 

2 



10 



-4 



/ 


':-^r.''" 


' ' ' z /d 
o' 


/: 




(IVzz')''^ ■■• 




- 








- 


'■ /^^ 


^L ~~~~~~~~^ 




; 


■ 






- 








V 


"1 , , , 




(b) 


\ 



0.5 



1 
m 



1.5 



Fig. 13. — Properties of the simple three-pole geometry used to estimate several factors in the MCC model, (a) Three-dimensional projection of selected field 
lines (gray curves), shown along with the two positive poles (filled circles) and the negative pole (open circle) on the surface, the sepai'ator field line (black solid 
curve), and outlines of the locations of the separatrix surfaces (dotted curves), (b) Plot that shows how the null-point height zo/d (dashed curve), the magnetic 
Jacobian factor (|CvaC;;|)''^ (dotted curve), and the constant Cl (solid curve) depend on the flux imbalance ratio m. The range of values for Cl used when 

fflffeMi(^8ffi¥li^#^'fi-(9M^§ W}*^*pSgi!fV^pWl^?'^^) a set of open field lines that originates from the right-hand positive pole, 
(3) a set of closed field lines that connects the left and center poles, and (4) a set of closed field lines that connects the center and 
right poles. There are two separatrix surfaces that delineate the boundaries between these domains: a vertical surface that spans 
the y-z plane and is defined by the condition x = 0, and the upper half of a prolate spheroidal surface centered on the origin. The 
separator field line is the intersection of the two separatrix surfaces, and for this model it is a semicircle in the y-z plane. 

In order to solve Equation ( |A3t we need to evaluate the exact position of the separator First, we locate its maximum height 
Zo by looking for the height of the magnetic null point along the vertical line denoted by x = and y = 0. We use Equations 
( IA4l )-( lA6b to solve for the magnitude of the magnetic field strength, but we do not worry about its absolute normalization. Along 
the vertical line in question, Bx = By = 0, and we find that 



B- ex 



(£/2 + z2)3/2 



m 



(A7) 



We set B^ = and search for a nontrivial solution for zo > 0. This is a cubic polynomial equation, and Figure [T3lb) shows the 
numerical solution for the ratio zo/d as a function of m. Solutions exist only for m < 2. Due to the symmetry in our assumed 
system, the separator field line is confined to the plane x = 0, and it subtends a semicircular shape for y ^ 0. Thus, the separator 



obeys y +z = Zq, its length is L = ttzo, and we can use the geometrical factor 5 = 1 in Equation 

We estimate the average value of B'^ along the separator by just computing its v alue at the maximum height {x = y- 
At this point, the field's parallel direction points along the y axis, so Equation ( IA3b can be written 



0, 



Bl = 



dBx dBr dB- dBx 

dx dz dx dz 



The cross-derivatives in the second term are found to be zero, and it can be shown that 



where 



B\ 



C„ = 



nd^ ^l^-C-l ' 

x — 2 m 

{x+ 1)5/2 " ^^ 



z = zo). 
(A8) 

(A9) 
(AlO) 
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r = 



\ — 2x m 



and X = (zo/d)^. The two dimensionless factors given in Equations (lAlOl ) and (lAlll i are related to Equation ( IA8b via 



dx 






C,, 



dB, 






a 



Figure [T3l b) shows the dimensionless factor (iCxiQ;!)'''^ as a function of the flux imbalance ratio m. 
The above model gives us the ability to write the average power dissipated (Equation dAll) ) as 



P = ^lCl 



*+ 



d^ 



dt 



where the dimensionless factors dependent on m have been collected into a single constant 



c.-^^VWM(^j) 



(All) 



(A12) 



(A13) 



(A14) 



Figure [T3l b) shows that Cl varies less strongly as a function of m than either of its components. In our models, we will not keep 
track of the individual m imbalance ratios for each reconnection event. Instead we adopt a range of values for Cl that spans the 
majority of the variation for many likely m values. The gray region in Figure [T3l b) shows this range of values; the lower limit is 
0.003, and the upper limit is the maximum valu e of .011. 

The other dimensionless constant in Equation (lAlb is 9]^. This parameter is a threshold ratio of the instantaneous current density 
to the characteristic current /*, and in the MCC model it is assumed that plasma instabiliti es (e.g., the ion -acoustic instability 
or tearing-mode instabilities) will limit the growth of the current to some fraction of /*. ILongco pe ( 1996) argued t hat 0l ^ 1 
was reasonable to expect , and h e ended up using 0l = 0.15 in the initial MCC models. However. iLongcope & Silval tl998) and 
ILongcope & Kankelborgid 19991) found that some situations appear to demand larger values of order ^l ~ 1- We will use the latter 
value, but we will keep in mind that the resulting heating rate may be an upper limit. 

To apply the heating rate derived above to our Monte Carlo models, we note that Equation ( fTSI l gives the ti me de rivative of 
magnetic flux that is being opened up in the simulation box, during each time step. In order to solve Equation (IA13l l, however, 
we also need to know the characteristic fluxes in the elements that are interacting, as well as their inter-element distances. Since 
many reconnection events may be occurring simultaneously in each time step, we must use averages taken over the box area. 
We also divide both sides of Equation ( IA13b by A in order to express the heating rate per unit area in terms of the variations in 
magnetic flux density. Thus, 



^lCi 



1l 

{d) 



dB 



dt 



where $1 = ^abs/N is the mean absolute flux per element in the simulation box, and 



^d) = \^, 




(A15) 



(A16) 



is the mean separation between flux elements. This is the form of the MCC energy flux used for the BONES results presented in 
Section[ 
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